{"cells":[{"cell_type":"markdown","metadata":{},"source":["# Iterated entanglement swapping\n","\n","**Download this notebook - {nb-download}`entanglement_swapping.ipynb`**"]},{"cell_type":"markdown","metadata":{},"source":["In this tutorial, we will focus on:\n","- designing circuits with mid-circuit measurement and conditional gates;\n","- utilising noise models in supported simulators."]},{"cell_type":"markdown","metadata":{},"source":["This example assumes the reader is familiar with the Qubit Teleportation and Entanglement Swapping protocols, and basic models of noise in quantum devices.\n","\n","To run this example, you will need `pytket`, `pytket-qiskit`, and `plotly` (installed via `pip`). To view the graphs, you will need an intallation of `plotly-orca`.\n","\n","Current quantum hardware fits into the NISQ (Noisy, Intermediate-Scale Quantum) regime. This noise cannot realistically be combatted using conventional error correcting codes, because of the lack of available qubits, noise levels exceeding the code thresholds, and very few devices available that can perform measurements and corrections mid-circuit. Analysis of how quantum algorithms perform under noisy conditions is a very active research area, as is finding ways to cope with it. Here, we will look at how well we can perform the Entanglement Swapping protocol with different noise levels.\n","\n","The Entanglement Swapping protocol requires two parties to share Bell pairs with a third party, who applies the Qubit Teleportation protocol to generate a Bell pair between the two parties. The Qubit Teleportation step requires us to be able to measure some qubits and make subsequent corrections to the remaining qubits. There are only a handful of simulators and devices that currently support this, with others restricted to only measuring the qubits at the end of the circuit.\n","\n","The most popular circuit model with conditional gates at the moment is that provided by the OpenQASM language. This permits a very restricted model of classical logic, where we can apply a gate conditionally on the exact value of a classical register. There is no facility in the current spec for Boolean logic or classical operations to apply any function to the value prior to the equality check. For example, Qubit Teleportation can be performed by the following QASM:\n","`OPENQASM 2.0;`\n","`include \"qelib1.inc\";`\n","`qreg a[2];`\n","`qreg b[1];`\n","`creg c[2];`\n","`// Bell state between Alice and Bob`\n","`h a[1];`\n","`cx a[1],b[0];`\n","`// Bell measurement of Alice's qubits`\n","`cx a[0],a[1];`\n","`h a[0];`\n","`measure a[0] -> c[0];`\n","`measure a[1] -> c[1];`\n","`// Correction of Bob's qubit`\n","`if(c==1) z b[0];`\n","`if(c==3) z b[0];`\n","`if(c==2) x b[0];`\n","`if(c==3) x b[0];`\n","\n","This corresponds to the following `pytket` code:"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from pytket import Circuit"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["qtel = Circuit()\n","alice = qtel.add_q_register(\"a\", 2)\n","bob = qtel.add_q_register(\"b\", 1)\n","data = qtel.add_c_register(\"d\", 2)"]},{"cell_type":"markdown","metadata":{},"source":["Bell state between Alice and Bob:"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["qtel.H(alice[1])\n","qtel.CX(alice[1], bob[0])"]},{"cell_type":"markdown","metadata":{},"source":["Bell measurement of Alice's qubits:"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["qtel.CX(alice[0], alice[1])\n","qtel.H(alice[0])\n","qtel.Measure(alice[0], data[0])\n","qtel.Measure(alice[1], data[1])"]},{"cell_type":"markdown","metadata":{},"source":["Correction of Bob's qubit:"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["qtel.X(bob[0], condition_bits=[data[0], data[1]], condition_value=2)\n","qtel.X(bob[0], condition_bits=[data[0], data[1]], condition_value=3)\n","qtel.Z(bob[0], condition_bits=[data[0], data[1]], condition_value=1)\n","qtel.Z(bob[0], condition_bits=[data[0], data[1]], condition_value=3)"]},{"cell_type":"markdown","metadata":{},"source":["So to demonstrate the Entanglement Swapping protocol, we just need to run this on one side of a Bell pair."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["es = Circuit()\n","ava = es.add_q_register(\"a\", 1)\n","bella = es.add_q_register(\"b\", 2)\n","charlie = es.add_q_register(\"c\", 1)\n","data = es.add_c_register(\"d\", 2)"]},{"cell_type":"markdown","metadata":{},"source":["Bell state between Ava and Bella:"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["es.H(ava[0])\n","es.CX(ava[0], bella[0])"]},{"cell_type":"markdown","metadata":{},"source":["Teleport `bella[0]` to `charlie[0]`:"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["tel_to_c = qtel.copy()\n","tel_to_c.rename_units({alice[0]: bella[0], alice[1]: bella[1], bob[0]: charlie[0]})\n","es.append(tel_to_c)"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["print(es.get_commands())"]},{"cell_type":"markdown","metadata":{},"source":["Let's start by running a noiseless simulation of this to verify that what we get looks like a Bell pair."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from pytket.extensions.qiskit import AerBackend"]},{"cell_type":"markdown","metadata":{},"source":["Connect to a simulator:"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["backend = AerBackend()"]},{"cell_type":"markdown","metadata":{},"source":["Make a ZZ measurement of the Bell pair:"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["bell_test = es.copy()\n","bell_test.Measure(ava[0], data[0])\n","bell_test.Measure(charlie[0], data[1])"]},{"cell_type":"markdown","metadata":{},"source":["Run the experiment:"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["bell_test = backend.get_compiled_circuit(bell_test)\n","from pytket.circuit.display import render_circuit_jupyter"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["render_circuit_jupyter(bell_test)\n","handle = backend.process_circuit(bell_test, n_shots=2000)\n","counts = backend.get_result(handle).get_counts()"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["print(counts)"]},{"cell_type":"markdown","metadata":{},"source":["This is good, we have got roughly 50/50 measurement results of 00 and 11 under the ZZ operator. But there are many other states beyond the Bell state that also generate this distribution, so to gain more confidence in our claim about the state we should make more measurements that also characterise it, i.e. perform state tomography.\n","\n","Here, we will demonstrate a naive approach to tomography that makes 3^n measurement circuits for an n-qubit state. More elaborate methods also exist."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from pytket.pauli import Pauli, QubitPauliString\n","from pytket.utils import append_pauli_measurement, probs_from_counts\n","from itertools import product\n","from scipy.linalg import lstsq, eigh\n","import numpy as np"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def gen_tomography_circuits(state, qubits, bits):\n","    # Yields {X, Y, Z}^n measurements in lexicographical order\n","    # Only measures qubits, storing the result in bits\n","    # (since we don't care about the ancilla qubits)\n","    assert len(qubits) == len(bits)\n","    for paulis in product([Pauli.X, Pauli.Y, Pauli.Z], repeat=len(qubits)):\n","        circ = state.copy()\n","        for qb, b, p in zip(qubits, bits, paulis):\n","            if p == Pauli.X:\n","                circ.H(qb)\n","            elif p == Pauli.Y:\n","                circ.V(qb)\n","            circ.Measure(qb, b)\n","        yield circ"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def run_tomography_circuits(state, qubits, bits, backend):\n","    circs = list(gen_tomography_circuits(state, qubits, bits))\n","    # Compile and run each circuit\n","    circs = backend.get_compiled_circuits(circs)\n","    handles = backend.process_circuits(circs, n_shots=2000)\n","    # Get the observed measurement probabilities\n","    probs_list = []\n","    for result in backend.get_results(handles):\n","        counts = result.get_counts()\n","        probs = probs_from_counts(counts)\n","        probs_list.append(probs)\n","    return probs_list"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def fit_tomography_outcomes(probs_list, n_qbs):\n","    # Define the density matrices for the basis states\n","    basis = dict()\n","    basis[(Pauli.X, 0)] = np.asarray([[0.5, 0.5], [0.5, 0.5]])\n","    basis[(Pauli.X, 1)] = np.asarray([[0.5, -0.5], [-0.5, 0.5]])\n","    basis[(Pauli.Y, 0)] = np.asarray([[0.5, -0.5j], [0.5j, 0.5]])\n","    basis[(Pauli.Y, 1)] = np.asarray([[0.5, 0.5j], [-0.5j, 0.5]])\n","    basis[(Pauli.Z, 0)] = np.asarray([[1, 0], [0, 0]])\n","    basis[(Pauli.Z, 1)] = np.asarray([[0, 0], [0, 1]])\n","    dim = 2**n_qbs\n","    # Define vector all_probs as a concatenation of probability vectors for each measurement (2**n x 3**n, 1)\n","    # Define matrix all_ops mapping a (vectorised) density matrix to a vector of probabilities for each measurement\n","    # (2**n x 3**n, 2**n x 2**n)\n","    all_probs = []\n","    all_ops = []\n","    for paulis, probs in zip(\n","        product([Pauli.X, Pauli.Y, Pauli.Z], repeat=n_qbs), probs_list\n","    ):\n","        prob_vec = []\n","        meas_ops = []\n","        for outcome in product([0, 1], repeat=n_qbs):\n","            prob_vec.append(probs.get(outcome, 0))\n","            op = np.eye(1, dtype=complex)\n","            for p, o in zip(paulis, outcome):\n","                op = np.kron(op, basis[(p, o)])\n","            meas_ops.append(op.reshape(1, dim * dim).conj())\n","        all_probs.append(np.vstack(prob_vec))\n","        all_ops.append(np.vstack(meas_ops))\n","    # Solve for density matrix by minimising || all_ops * dm - all_probs ||\n","    dm, _, _, _ = lstsq(np.vstack(all_ops), np.vstack(all_probs))\n","    dm = dm.reshape(dim, dim)\n","    # Make density matrix positive semi-definite\n","    v, w = eigh(dm)\n","    for i in range(dim):\n","        if v[i] < 0:\n","            for j in range(i + 1, dim):\n","                v[j] += v[i] / (dim - (i + 1))\n","            v[i] = 0\n","    dm = np.zeros([dim, dim], dtype=complex)\n","    for j in range(dim):\n","        dm += v[j] * np.outer(w[:, j], np.conj(w[:, j]))\n","    # Normalise trace of density matrix\n","    dm /= np.trace(dm)\n","    return dm"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["probs_list = run_tomography_circuits(\n","    es, [ava[0], charlie[0]], [data[0], data[1]], backend\n",")\n","dm = fit_tomography_outcomes(probs_list, 2)\n","print(dm.round(3))"]},{"cell_type":"markdown","metadata":{},"source":["This is very close to the true density matrix for a pure Bell state. We can attribute the error here to the sampling error since we only take 2000 samples of each measurement circuit.\n","\n","To quantify exactly how similar it is to the correct density matrix, we can calculate the fidelity."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from scipy.linalg import sqrtm"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def fidelity(dm0, dm1):\n","    # Calculate the fidelity between two density matrices\n","    sq0 = sqrtm(dm0)\n","    sq1 = sqrtm(dm1)\n","    return np.linalg.norm(sq0.dot(sq1)) ** 2"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["bell_state = np.asarray(\n","    [\n","        [0.5, 0, 0, 0.5],\n","        [0, 0, 0, 0],\n","        [0, 0, 0, 0],\n","        [0.5, 0, 0, 0.5],\n","    ]\n",")\n","print(fidelity(dm, bell_state))"]},{"cell_type":"markdown","metadata":{},"source":["This high fidelity is unsurprising since we have a completely noiseless simulation. So the next step is to add some noise to the simulation and observe how the overall fidelity is affected. The `AerBackend` wraps around the Qiskit Aer simulator and can pass on any `qiskit.providers.aer.noise.NoiseModel` to the simulator. Let's start by adding some uniform depolarising noise to each CX gate and some uniform measurement error."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from qiskit_aer.noise import NoiseModel, depolarizing_error, ReadoutError"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def make_noise_model(dep_err_rate, ro_err_rate, qubits):\n","    # Define a noise model that applies uniformly to the given qubits\n","    model = NoiseModel()\n","    dep_err = depolarizing_error(dep_err_rate, 2)\n","    ro_err = ReadoutError(\n","        [[1 - ro_err_rate, ro_err_rate], [ro_err_rate, 1 - ro_err_rate]]\n","    )\n","    # Add depolarising error to CX gates between any qubits (implying full connectivity)\n","    for i, j in product(qubits, repeat=2):\n","        if i != j:\n","            model.add_quantum_error(dep_err, [\"cx\"], [i, j])\n","    # Add readout error for each qubit\n","    for i in qubits:\n","        model.add_readout_error(ro_err, qubits=[i])\n","    return model"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["test_model = make_noise_model(0.03, 0.05, range(4))\n","backend = AerBackend(noise_model=test_model)\n","probs_list = run_tomography_circuits(\n","    es, [ava[0], charlie[0]], [data[0], data[1]], backend\n",")\n","dm = fit_tomography_outcomes(probs_list, 2)\n","print(dm.round(3))\n","print(fidelity(dm, bell_state))"]},{"cell_type":"markdown","metadata":{},"source":["Despite the very small circuit and the relatively small error rates, the fidelity of the final state has reduced considerably.\n","\n","As far as circuits go, the entanglement swapping protocol is little more than a toy example and is nothing close to the scale of circuits for most interesting quantum computational problems. However, it is possible to iterate the protocol many times to build up a larger computation, allowing us to see the impact of the noise at different scales."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from pytket import OpType\n","from plotly.graph_objects import Scatter, Figure"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def iterated_entanglement_swap(n_iter):\n","    # Iterate the entanglement swapping protocol n_iter times\n","    it_es = Circuit()\n","    ava = it_es.add_q_register(\"a\", 1)\n","    bella = it_es.add_q_register(\"b\", 2)\n","    charlie = it_es.add_q_register(\"c\", 1)\n","    data = it_es.add_c_register(\"d\", 2)\n","\n","    # Start with an initial Bell state\n","    it_es.H(ava[0])\n","    it_es.CX(ava[0], bella[0])\n","    for i in range(n_iter):\n","        if i % 2 == 0:\n","            # Teleport bella[0] to charlie[0] to give a Bell pair between ava[0] and charlier[0]\n","            tel_to_c = qtel.copy()\n","            tel_to_c.rename_units(\n","                {alice[0]: bella[0], alice[1]: bella[1], bob[0]: charlie[0]}\n","            )\n","            it_es.append(tel_to_c)\n","            it_es.add_gate(OpType.Reset, [bella[0]])\n","            it_es.add_gate(OpType.Reset, [bella[1]])\n","        else:\n","            # Teleport charlie[0] to bella[0] to give a Bell pair between ava[0] and bella[0]\n","            tel_to_b = qtel.copy()\n","            tel_to_b.rename_units(\n","                {alice[0]: charlie[0], alice[1]: bella[1], bob[0]: bella[0]}\n","            )\n","            it_es.append(tel_to_b)\n","            it_es.add_gate(OpType.Reset, [bella[1]])\n","            it_es.add_gate(OpType.Reset, [charlie[0]])\n","    # Return the circuit and the qubits expected to share a Bell pair\n","    if n_iter % 2 == 0:\n","        return it_es, [ava[0], bella[0]]\n","    else:\n","        return it_es, [ava[0], charlie[0]]"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def iterated_noisy_experiment(dep_err_rate, ro_err_rate, max_iter):\n","    # Set up the noisy simulator with the given error rates\n","    test_model = make_noise_model(dep_err_rate, ro_err_rate, range(4))\n","    backend = AerBackend(noise_model=test_model)\n","    # Estimate the fidelity after n iterations, from 0 to max_iter (inclusive)\n","    fid_list = []\n","    for i in range(max_iter + 1):\n","        it_es, qubits = iterated_entanglement_swap(i)\n","        probs_list = run_tomography_circuits(it_es, qubits, [data[0], data[1]], backend)\n","        dm = fit_tomography_outcomes(probs_list, 2)\n","        fid = fidelity(dm, bell_state)\n","        fid_list.append(fid)\n","    return fid_list"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["fig = Figure()\n","fig.update_layout(\n","    title=\"Iterated Entanglement Swapping under Noise (dep_err = 0.03)\",\n","    xaxis_title=\"Iterations\",\n","    xaxis=dict(range=[0, 10]),\n","    yaxis_title=\"Fidelity\",\n",")\n","iter_range = np.arange(11)\n","for i in range(7):\n","    fids = iterated_noisy_experiment(0.03, 0.025 * i, 10)\n","    plot_data = Scatter(\n","        x=iter_range, y=fids, name=\"ro_err=\" + str(np.round(0.025 * i, 3))\n","    )\n","    fig.add_trace(plot_data)\n","try:\n","    fig.show(renderer=\"svg\")\n","except ValueError as e:\n","    print(e)  # requires plotly-orca"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["fig = Figure()\n","fig.update_layout(\n","    title=\"Iterated Entanglement Swapping under Noise (ro_err = 0.05)\",\n","    xaxis_title=\"Iterations\",\n","    xaxis=dict(range=[0, 10]),\n","    yaxis_title=\"Fidelity\",\n",")\n","iter_range = np.arange(11)\n","for i in range(9):\n","    fids = iterated_noisy_experiment(0.01 * i, 0.05, 10)\n","    plot_data = Scatter(\n","        x=iter_range, y=fids, name=\"dep_err=\" + str(np.round(0.01 * i, 3))\n","    )\n","    fig.add_trace(plot_data)\n","try:\n","    fig.show(renderer=\"svg\")\n","except ValueError as e:\n","    print(e)  # requires plotly-orca"]},{"cell_type":"markdown","metadata":{},"source":["These graphs are not very surprising, but are still important for seeing that the current error rates of typical NISQ devices become crippling for fidelities very quickly after repeated mid-circuit measurements and corrections (even with this overly-simplified model with uniform noise and no crosstalk or higher error modes). This provides good motivation for the adoption of error mitigation techniques, and for the development of new techniques that are robust to errors in mid-circuit measurements."]},{"cell_type":"markdown","metadata":{},"source":["Exercises:\n","- Vary the fixed noise levels to compare how impactful the depolarising and measurement errors are.\n","- Add extra noise characteristics to the noise model to obtain something that more resembles a real device. Possible options include adding error during the reset operations, extending the errors to be non-local, or constructing the noise model from a device's calibration data.\n","- Change the circuit from iterated entanglement swapping to iterated applications of a correction circuit from a simple error-correcting code. Do you expect this to be more sensitive to depolarising errors from unitary gates or measurement errors?"]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.4"}},"nbformat":4,"nbformat_minor":2}