{"cells":[{"cell_type":"markdown","metadata":{},"source":["# Expectation values\n","\n","**Download this notebook - {nb-download}`expectation_value_example.ipynb`**"]},{"cell_type":"markdown","metadata":{},"source":["Given a circuit generating a quantum state $\\lvert \\psi \\rangle$, it is very common to have an operator $H$ and ask for the expectation value $\\langle \\psi \\vert H \\vert \\psi \\rangle$. A notable example is in quantum computational chemistry, where $\\lvert \\psi \\rangle$ encodes the wavefunction for the electronic state of a small molecule, and the energy of the molecule can be derived from the expectation value with respect to the molecule's Hamiltonian operator $H$.<br>\n","<br>\n","This example uses this chemistry scenario to demonstrate the overall procedure for using `pytket` to perform advanced high-level procedures. We build on top of topics covered by several other example notebooks, including circuit generation, optimisation, and using different backends.<br>\n","<br>\n","There is limited built-in functionality in `pytket` for obtaining expectation values from circuits. This is designed to encourage users to consider their needs for parallelising the processing of circuits, manipulating results (e.g. filtering, adjusting counts to mitigate errors, and other forms of data processing), or more advanced schemes for grouping the terms of the operator into measurement circuits. For this example, suppose that we want to focus on reducing the queueing time for IBM device backends, and filter our shots to eliminate some detected errors.<br>\n","<br>\n","This notebook makes use of the Qiskit and ProjectQ backend modules `pytket_qiskit` and `pytket_projectq`, as well as the electronic structure module `openfermion`, all three of which should first be installed via `pip`.<br>\n","<br>\n","We will start by generating an ansatz and Hamiltonian for the chemical of interest. Here, we are just using a simple model of $\\mathrm{H}_2$ with four qubits representing the occupation of four spin orbitals."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from pytket import Circuit, Qubit, Bit\n","from sympy import symbols"]},{"cell_type":"markdown","metadata":{},"source":["Generate ansatz and Hamiltonian:"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["ansatz = Circuit()\n","qubits = ansatz.add_q_register(\"q\", 4)\n","args = symbols(\"a0 a1 a2 a3 a4 a5 a6 a7\")\n","for i in range(4):\n","    ansatz.Ry(args[i], qubits[i])\n","for i in range(3):\n","    ansatz.CX(qubits[i], qubits[i + 1])\n","for i in range(4):\n","    ansatz.Ry(args[4 + i], qubits[i])\n","ansatz.measure_all()"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["for command in ansatz:\n","    print(command)"]},{"cell_type":"markdown","metadata":{},"source":["In reality, you would use an expectation value calculation as the objective function for a classical optimisation routine to determine the parameter values for the ground state. For the purposes of this notebook, we will use some predetermined values for the ansatz, already optimised for $\\mathrm{H}_2$."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["arg_values = [\n","    7.17996183e-02,\n","    2.95442468e-08,\n","    1.00000015e00,\n","    1.00000086e00,\n","    9.99999826e-01,\n","    1.00000002e00,\n","    9.99999954e-01,\n","    1.13489747e-06,\n","]"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["ansatz.symbol_substitution(dict(zip(args, arg_values)))"]},{"cell_type":"markdown","metadata":{},"source":["We can use for example the openfermion library to express an Hamiltonian as a sum of tensors of paulis."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["import openfermion as of"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["hamiltonian = (\n","    -0.0970662681676282 * of.QubitOperator(\"\")\n","    + -0.045302615503799284 * of.QubitOperator(\"X0 X1 Y2 Y3\")\n","    + 0.045302615503799284 * of.QubitOperator(\"X0 Y1 Y2 X3\")\n","    + 0.045302615503799284 * of.QubitOperator(\"Y0 X1 X2 Y3\")\n","    + -0.045302615503799284 * of.QubitOperator(\"Y0 Y1 X2 X3\")\n","    + 0.17141282644776884 * of.QubitOperator(\"Z0\")\n","    + 0.16868898170361213 * of.QubitOperator(\"Z0 Z1\")\n","    + 0.12062523483390425 * of.QubitOperator(\"Z0 Z2\")\n","    + 0.16592785033770352 * of.QubitOperator(\"Z0 Z3\")\n","    + 0.17141282644776884 * of.QubitOperator(\"Z1\")\n","    + 0.16592785033770352 * of.QubitOperator(\"Z1 Z2\")\n","    + 0.12062523483390425 * of.QubitOperator(\"Z1 Z3\")\n","    + -0.22343153690813597 * of.QubitOperator(\"Z2\")\n","    + 0.17441287612261608 * of.QubitOperator(\"Z2 Z3\")\n","    + -0.22343153690813597 * of.QubitOperator(\"Z3\")\n",")"]},{"cell_type":"markdown","metadata":{},"source":["This can be converted into pytket's QubitPauliOperator type.<br>\n","<br>\n","The OpenFermion `QubitOperator` class represents the operator by its decomposition into a linear combination of Pauli operators (tensor products of the $I$, $X$, $Y$, and $Z$ matrices).<br>\n","<br>\n","A `QubitPauliString` is a sparse representation of a Pauli operator with support over some subset of qubits."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from pytket.pauli import Pauli, QubitPauliString\n","from pytket.utils.operators import QubitPauliOperator"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["pauli_sym = {\"I\": Pauli.I, \"X\": Pauli.X, \"Y\": Pauli.Y, \"Z\": Pauli.Z}"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def qps_from_openfermion(paulis):\n","    \"\"\"Convert OpenFermion tensor of Paulis to pytket QubitPauliString.\"\"\"\n","    qlist = []\n","    plist = []\n","    for q, p in paulis:\n","        qlist.append(Qubit(q))\n","        plist.append(pauli_sym[p])\n","    return QubitPauliString(qlist, plist)"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def qpo_from_openfermion(openf_op):\n","    \"\"\"Convert OpenFermion QubitOperator to pytket QubitPauliOperator.\"\"\"\n","    tk_op = dict()\n","    for term, coeff in openf_op.terms.items():\n","        string = qps_from_openfermion(term)\n","        tk_op[string] = coeff\n","    return QubitPauliOperator(tk_op)"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["hamiltonian_op = qpo_from_openfermion(hamiltonian)"]},{"cell_type":"markdown","metadata":{},"source":["We can simulate this exactly using a statevector simulator like ProjectQ. This has a built-in method for fast calculations of expectation values that works well for small examples like this."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from pytket.extensions.projectq import ProjectQBackend"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["backend = ProjectQBackend()\n","ideal_energy = backend.get_operator_expectation_value(ansatz, hamiltonian_op)\n","print(ideal_energy)"]},{"cell_type":"markdown","metadata":{},"source":["Ideally the state generated by this ansatz will only span the computational basis states with exactly two of the four qubits in state $\\lvert 1 \\rangle$. This is because these basis states correspond to two electrons being present in the molecule.<br>\n","<br>\n","This ansatz is a hardware-efficient model that is designed to explore a large portion of the Hilbert space with relatively few entangling gates. Unfortunately, with this much freedom, it will regularly generate states that have no physical interpretation such as states spanning multiple basis states corresponding to different numbers of electrons in the system (which we assume is fixed and conserved).<br>\n","<br>\n","We can mitigate this by using a syndrome qubit that calculates the parity of the other qubits. Post-selecting this syndrome with $\\langle 0 \\rvert$ will project the remaining state onto the subspace of basis states with even parity, increasing the likelihood the observed state will be a physically admissible state.<br>\n","<br>\n","Even if the ansatz parameters are tuned to give a physical state, real devices have noise and imperfect gates, so in practice we may also measure bad states with a small probability. If this syndrome qubit is measured as 1, it means an error has definitely occurred, so we should discard the shot."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["syn = Qubit(\"synq\", 0)\n","syn_res = Bit(\"synres\", 0)\n","ansatz.add_qubit(syn)\n","ansatz.add_bit(syn_res)\n","for qb in qubits:\n","    ansatz.CX(qb, syn)\n","ansatz.Measure(syn, syn_res)"]},{"cell_type":"markdown","metadata":{},"source":["Using this, we can define a filter function which removes the shots which the syndrome qubit detected as erroneous. `BackendResult` objects allow retrieval of shots in any bit order, so we can retrieve the `synres` results separately and use them to filter the shots from the remaining bits. The Backends example notebook describes this in more detail."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from collections import Counter"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def filter_shots(backend_result, syn_res_bit):\n","    bits = sorted(backend_result.get_bitlist())\n","    bits.remove(syn_res_bit)\n","    syn_shots = backend_result.get_shots([syn_res])[:, 0]\n","    main_shots = backend_result.get_shots(bits)\n","    return main_shots[syn_shots == 0]"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def filter_counts(backend_result, syn_res_bit):\n","    bits = sorted(backend_result.get_bitlist())\n","    syn_index = bits.index(syn_res_bit)\n","    counts = backend_result.get_counts()\n","    filtered_counts = Counter()\n","    for readout, count in counts.items():\n","        if readout[syn_index] == 0:\n","            filtered_readout = tuple(v for i, v in enumerate(readout) if i != syn_index)\n","            filtered_counts[filtered_readout] += count\n","    return filtered_counts"]},{"cell_type":"markdown","metadata":{},"source":["Depending on which backend we will be using, we will need to compile each circuit we run to conform to the gate set and connectivity constraints. We can define a compilation pass for each backend that optimises the circuit and maps it onto the backend's gate set and connectivity constraints. We don't expect this to change our circuit too much as it is already near-optimal."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from pytket.passes import OptimisePhaseGadgets, SequencePass"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def compiler_pass(backend):\n","    return SequencePass([OptimisePhaseGadgets(), backend.default_compilation_pass()])"]},{"cell_type":"markdown","metadata":{},"source":["Given the full statevector, the expectation value can be calculated simply by matrix multiplication. However, with a real quantum system, we cannot observe the full statevector directly. Fortunately, the Pauli decomposition of the operator gives us a sequence of measurements we should apply to obtain the relevant information to reconstruct the expectation value.<br>\n","<br>\n","The utility method `append_pauli_measurement` takes a single term of a `QubitPauliOperator` (a `QubitPauliString`) and appends measurements in the corresponding bases to obtain the expectation value for that particular Pauli operator. We will want to make a new `Circuit` object for each of the measurements we wish to observe.<br>\n"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from pytket.predicates import CompilationUnit\n","from pytket.utils import append_pauli_measurement"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def gen_pauli_measurement_circuits(state_circuit, compiler_pass, operator):\n","    # compile main circuit once\n","    state_cu = CompilationUnit(state_circuit)\n","    compiler_pass.apply(state_cu)\n","    compiled_state = state_cu.circuit\n","    final_map = state_cu.final_map\n","    # make a measurement circuit for each pauli\n","    pauli_circuits = []\n","    coeffs = []\n","    energy = 0\n","    for p, c in operator.terms.items():\n","        if p == ():\n","            # constant term\n","            energy += c\n","        else:\n","            # make measurement circuits and compile them\n","            pauli_circ = Circuit(state_circuit.n_qubits - 1)  # ignore syndrome qubit\n","            append_pauli_measurement(qps_from_openfermion(p), pauli_circ)\n","            pauli_cu = CompilationUnit(pauli_circ)\n","            compiler_pass.apply(pauli_cu)\n","            pauli_circ = pauli_cu.circuit\n","            init_map = pauli_cu.initial_map\n","            # map measurements onto the placed qubits from the state\n","            rename_map = {\n","                i: final_map[o] for o, i in init_map.items() if o in final_map\n","            }\n","            pauli_circ.rename_units(rename_map)\n","            state_and_measure = compiled_state.copy()\n","            state_and_measure.append(pauli_circ)\n","            pauli_circuits.append(state_and_measure)\n","            coeffs.append(c)\n","    return pauli_circuits, coeffs, energy"]},{"cell_type":"markdown","metadata":{},"source":["We can now start composing these together to get our generalisable expectation value function. Passing all of our circuits to `process_circuits` allows them to be submitted to IBM Quantum devices at the same time, giving substantial savings in overall queueing time. Since the backend will cache any results from `Backend.process_circuits`, we will remove the results when we are done with them to prevent memory bloating when this method is called many times."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from pytket.utils import expectation_from_shots, expectation_from_counts"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["def expectation_value(state_circuit, operator, backend, n_shots):\n","    if backend.supports_expectation:\n","        circuit = state_circuit.copy()\n","        compiled_circuit = backend.get_compiled_circuit(circuit)\n","        return backend.get_operator_expectation_value(\n","            compiled_circuit, qpo_from_openfermion(operator)\n","        )\n","    elif backend.supports_shots:\n","        syn_res_index = state_circuit.bit_readout[syn_res]\n","        pauli_circuits, coeffs, energy = gen_pauli_measurement_circuits(\n","            state_circuit, compiler_pass(backend), operator\n","        )\n","        handles = backend.process_circuits(pauli_circuits, n_shots=n_shots)\n","        for handle, coeff in zip(handles, coeffs):\n","            res = backend.get_result(handle)\n","            filtered = filter_shots(res, syn_res)\n","            energy += coeff * expectation_from_shots(filtered)\n","            backend.pop_result(handle)\n","        return energy\n","    elif backend.supports_counts:\n","        syn_res_index = state_circuit.bit_readout[syn_res]\n","        pauli_circuits, coeffs, energy = gen_pauli_measurement_circuits(\n","            state_circuit, compiler_pass(backend), operator\n","        )\n","        handles = backend.process_circuits(pauli_circuits, n_shots=n_shots)\n","        for handle, coeff in zip(handles, coeffs):\n","            res = backend.get_result(handle)\n","            filtered = filter_counts(res, syn_res)\n","            energy += coeff * expectation_from_counts(filtered)\n","            backend.pop_result(handle)\n","        return energy\n","    else:\n","        raise NotImplementedError(\"Implementation for state to be written\")"]},{"cell_type":"markdown","metadata":{},"source":["...and then run it for our ansatz. `AerBackend` supports faster expectation value from snapshopts (using the `AerBackend.get_operator_expectation_value` method), but this only works when all the qubits in the circuit are default register qubits that go up from 0. So we will need to rename `synq`."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["from pytket.extensions.qiskit import IBMQEmulatorBackend, AerBackend"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["ansatz.rename_units({Qubit(\"synq\", 0): Qubit(\"q\", 4)})"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["print(expectation_value(ansatz, hamiltonian, AerBackend(), 8000))\n","# Try replacing IBMQEmulatorBackend with IBMQBackend to submit the circuits to a real IBM Quantum device.\n","print(expectation_value(ansatz, hamiltonian, IBMQEmulatorBackend(\"ibmq_manila\"), 8000))"]},{"cell_type":"markdown","metadata":{},"source":["For basic practice with using pytket backends and their results, try editing the code here to:<br>\n","* Extend `expectation_value` to work with statevector backends (e.g. `AerStateBackend`)<br>\n","* Remove the row filtering from `filter_shots` and see the effect on the expectation value on a noisy simulation/device<br>\n","* Adapt `filter_shots` to be able to filter a counts dictionary and adapt `expectation_value` to calulate the result using the counts summary from the backend (`pytket.utils.expectation_from_counts` will be useful here)"]}],"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}