jupyter_notebooks/Tutorials/algorithms/ErrorGeneratorPropagation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a2bf624d-6098-46e5-820e-e775e4fb41f0",
   "metadata": {},
   "source": [
    "# Error Generator Propagation\n",
    "In this tutorial we will provide an overview of the core functionality available through pyGSTi's error generator propagation module.\n",
    "\n",
    "Error generator propagation is a technique which leverages the analytical properties of the error generator formalism to enable efficient forward simulation by propagating general markovian error channels through Clifford circuits. Efficiency of this technique relies on two conditions:\n",
    "\n",
    "- Sparsity: At most a polynomial number of error generator rates (in the number of qubits) can be nonzero for any given circuit layer.\n",
    "- Clifford-only: The propagation of error generators relies on the analytic properties of the elementary error generators when conjugated by cliffords.\n",
    "\n",
    "That is pretty much it though. Coherent errors, non-unital errors (e.g. amplitude damping), dephasing, all fair game. Practically there is a third requirement as well and that is that the error generator rates are relatively small. The larger the error generator rates, the higher-order the approximation you'll require (BCH and/or taylor series) to achieve a given precision target when using the functionality described herein for efficiently performing strong simulation in the error generator propagation framework. \n",
    "\n",
    "Please note: The implementation of the error generator propagation framework in pyGSTi requires the `stim` python package, so please ensure this is installed in your environment before proceeding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "472bba68-9e69-4379-aa3c-7e062128ef7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pygsti\n",
    "import stim\n",
    "from pygsti.tools import errgenproptools as eprop\n",
    "from pygsti.tools.lindbladtools import random_error_generator_rates\n",
    "from pygsti.errorgenpropagation.errorpropagator import ErrorGeneratorPropagator\n",
    "from pygsti.errorgenpropagation.localstimerrorgen import LocalStimErrorgenLabel as _LSE"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a4682fa-f5e3-4a44-a596-e095f2cc4890",
   "metadata": {},
   "source": [
    "To begin we need an error model, and particularly one parameterized using error generators (or otherwise capable of outputing error generators for a circuit layer). For this tutorial we'll work with a 4-qubit crosstalk-free model for a gate set consisting of $\\pi/2$ rotations about X and Y on each qubit, and a two-qubit CPHASE gate. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae464763-b211-4bf4-a85c-47465aa7187b",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_qubits = 4\n",
    "gate_names = ['Gcphase', 'Gxpi2', 'Gypi2']\n",
    "availability = {'Gcphase':[(0,1), (1,2), (2,3), (3,0)]}\n",
    "pspec = pygsti.processors.QubitProcessorSpec(num_qubits, gate_names, availability=availability)\n",
    "target_model = pygsti.models.create_crosstalk_free_model(processor_spec = pspec)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f24def6b-b268-4a62-bbe5-37d0c188c15b",
   "metadata": {},
   "source": [
    "Now that we have a target model we'll also want a noisy model to simulate as well. For this example we'll randomly sample a weight-2 H+S (coherent + pauli stochastic) error model, but the error generator propagation framework can also handle C and A error generators as well (i.e. general lindbladian errors). \n",
    "The specific specification we'll need for the model construction routine we're about to use is a dictionary whose keys are gate labels. Each value of this dictionary is itself a dictionary whose keys are elementary error generator labels, and whose values are error generator rates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e205a715-c231-40a1-82bf-27cf303e8820",
   "metadata": {},
   "outputs": [],
   "source": [
    "qubit_labels = pspec.qubit_labels\n",
    "error_rates_dict = {}\n",
    "for gate, availability in pspec.availability.items():\n",
    "    n = pspec.gate_num_qubits(gate)\n",
    "    if availability == 'all-edges':\n",
    "        assert(n == 1), \"Currently require all 2-qubit gates have a specified availability!\"\n",
    "        qubits_for_gate = qubit_labels\n",
    "    else:\n",
    "        qubits_for_gate = availability  \n",
    "    for qs in qubits_for_gate:\n",
    "        label = pygsti.baseobjs.Label(gate, qs)\n",
    "        # Sample error rates.\n",
    "        error_rates_dict[label] = random_error_generator_rates(num_qubits=n, errorgen_types=('H', 'S'), label_type='local', seed=1234)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ebacf47-54ba-4cd1-b14c-26f1d0516d29",
   "metadata": {},
   "outputs": [],
   "source": [
    "error_model = pygsti.models.create_crosstalk_free_model(pspec, lindblad_error_coeffs=error_rates_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab4dc617-6d7c-410c-a976-282a169c8bdf",
   "metadata": {},
   "source": [
    "We'll also need an example circuit for the rest of our examples, so will construct one at random."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e88638b-cc35-49be-8976-ee658d3201a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = pygsti.algorithms.randomcircuit.create_random_circuit(pspec, 3, sampler='edgegrab', samplerargs=[0.4,], rand_state=12345)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "145d7875-599a-447f-b2ac-b2fc702d6dda",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "192d74b6-3f36-499f-9db6-335a01b87c3f",
   "metadata": {},
   "source": [
    "## Basic Propagation\n",
    "In this section we'll introduce the basic syntax of the `ErrorGeneratorPropagator` class and usage of the basic error generator propagation functionality.\n",
    "Our first step will be to create an instance of the `ErrorGeneratorPropagator` class. This is as simple as passing in our error model into the constructor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c66532cd-876e-4002-89d2-b47eaafb5bf0",
   "metadata": {},
   "outputs": [],
   "source": [
    "errorgen_propagator = ErrorGeneratorPropagator(error_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fc2f197-6351-4c16-b351-7e539f1db839",
   "metadata": {},
   "source": [
    "The very first thing we can do is propagate the error generators for each circuit layer to the end of the circuit. This is done using the `propagate_errorgens` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90dc2734-2d17-4555-b32c-ac307967c21e",
   "metadata": {},
   "outputs": [],
   "source": [
    "propagated_errorgen_layers = errorgen_propagator.propagate_errorgens(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "044babe6-675d-4145-bb68-bc3e5e80efff",
   "metadata": {},
   "source": [
    "The output of this method is a list of dictionaries, one for each original error generator layer in the circuit, containing an updated set of elementary error generator coefficients and rates corresponding to the result of propagating each error generator through the circuit. Note this list is returned in circuit ordering, so there is a one-to-one correspondence between the position an error generator appears in the original circuit and where it appears in this final list.\n",
    "\n",
    "So, to see the result of propagating the error generator corresponding to the noise induced after the first layer of gates to the very end we could query this list as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b43a4a84-2795-440a-9db4-9d80cfceea6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(propagated_errorgen_layers[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "981a3cb7-a629-44cd-bb17-bc52dd74e9a2",
   "metadata": {},
   "source": [
    "There are a few things worth noting at this point. First, we stated we'd be looking at the output of propagating the *first* circuit layer to the end but we indexed into the *second* position of the final list, what gives? This is because by default the `propagate_errorgens` method prepends and appends the error generator layers corresponding to state preparation and measurement respectively *before* beginning the propagation. As such the first layer in the final output corresponds to the error generator associated with state prep, and the final one with measurement. We never actually specified error generator rates for the SPAM, so you'll notice the corresponding dictionaries in the final output are both empty in this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7710d39-9574-41aa-a922-c95f7bb8225f",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(propagated_errorgen_layers[0])\n",
    "print(propagated_errorgen_layers[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85b38f4f-9797-416b-b04d-395caed0438b",
   "metadata": {},
   "source": [
    "To change this behavior so that the SPAM layers are not included you can set the optional kwarg `include_spam` to `FALSE` in `propgate_errorgens` and other related methods."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1076dfa4-da39-4c3d-9466-0de50c415521",
   "metadata": {},
   "source": [
    "The next things worth noting are the keys of the final dictionary. Notice that the basis element labels for each of the elementary error generator coefficient labels are instances of `stim.PauliString`, very much unlike the other elementary error generator labels used in pyGSTi. These labels are instances of the class `LocalStimErrorgenLabel`, a specialized label class with additional metadata and methods used throughout the error generator propagation framework. For applications where you need to take the output of this module and utilize it elsewhere in pyGSTi you can utilize the `to_local_eel` and `to_global_eel` methods of the `LocalStimErrorgenLabel` class to convert these into instances of `LocalElementaryErrorgenLabel` and `GlobalElementaryErrorgenLabel`, respectively, for use within other parts of pyGSTi."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af5abf19-7c96-4c06-aa41-4a2b47890205",
   "metadata": {},
   "source": [
    "While the output of `propgate_errorgens` is in and of itself incredibly useful, often we want to know more about how specific errors have been transformed by propagation through the circuit. Fortunately the analytic structure of error generator propagation through a clifford operation is such that it acts as a generalized permutation of each elementary error generator within it's sector (i.e. propagation can't in and of itself map H errors to anything other than H errors, for example). To view the input-output corresponding to the transformation of each error generator we can use the `errorgen_transform_map` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4b3ab03-399f-43dd-99d0-957bbf219ad5",
   "metadata": {},
   "outputs": [],
   "source": [
    "errorgen_transform_map = errorgen_propagator.errorgen_transform_map(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eeddc38f-78cc-46e9-ae36-c506a352c7a5",
   "metadata": {},
   "source": [
    "This method returns a dictionary with the following structure: Keys are tuples of the form (<original_errorgen_label>, <layer_index>), and values are of the form (<final_errorgen_label>, <overall_phase>), where overall_phase corresponds to the overall sign accumulated on the final error generator rate as a result of propagation. So, for example, we can see that as a result of propagation through the circuit the H(XIII) error generator at circuit layer 1 is mapped to an H(ZIII) error generator accruing and overall phase of -1.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5600efa-93a4-4100-93f7-43173c83f948",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(errorgen_transform_map[(_LSE('H', [stim.PauliString('XIII')]), 1)])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16b35551-56e1-4cb9-ab69-9e0bb71205d6",
   "metadata": {},
   "source": [
    "For some purposes it can be useful to go another step further and identity which gate a particular error might be associated with in the original error model. For this purpose `ErrorGeneratorPropagator` has a helper method available called `errorgen_gate_contributors`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "432a7f5d-d27b-4f12-a262-810609318cf3",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(errorgen_propagator.errorgen_gate_contributors(_LSE('H', [stim.PauliString('XIII')]), c, layer_idx=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd1e42b3-2309-4c67-9d78-67ac5510bbb9",
   "metadata": {},
   "source": [
    "Here this method returns the fact that in our particular error model the only gate at layer index 1 which could have contributed this particular error generator was the 'Gxpi2' gate acting on qubit 0. In some error models it may be possible for multiple gates to contribute to a particular rate, in which case this method should return all such gates."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90c2be36-4333-40b7-b3a4-742bc60623b1",
   "metadata": {},
   "source": [
    "## BCH Approximation\n",
    "In the previous section we showed how to use the `ErrorGeneratorPropagator` class to transform a circuit with a series of post-gate error generators into an equivalent representation of this noisy circuit with instead a series of post-circuit error generator layers. What if we want a single effective end-of-circuit error generator which approximates the overall action of the composition of each of the propagated error generators? To do so the `ErrorGeneratorPropagator` class supports the option to iteratively apply the BCH approximation at various orders to perform this recombination.\n",
    "\n",
    "The main method for performing propagation together with the iterative application of the BCH approximation is called `propagate_errorgens_bch`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ee0db5d-972c-4367-96da-f390fe1ac56e",
   "metadata": {},
   "outputs": [],
   "source": [
    "propagated_errorgen_layer_first_order = errorgen_propagator.propagate_errorgens_bch(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c5f5c55-fbc1-4371-8f7e-d274c6c361bc",
   "metadata": {},
   "source": [
    "As before this method propagated all of a circuits error generator layers to the very end, but follows this up with an iterative application of the BCH approximation resulting in a single final error generator. Without any additional optional arguments specified this uses the first-order BCH approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0d58c38-f3ce-494d-8daa-cbb5ef8d66b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(propagated_errorgen_layer_first_order)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7b85f44-fe5e-4d2e-9555-8cae50e8f03b",
   "metadata": {},
   "source": [
    "This method supports a number of additional arguments beyond those already for `propagate_errorgens`:\n",
    "- `bch_order`: An integer from 1 to 5 specifying the order of the BCH approximation to apply (5 is the current maximum). Note that the computational cost of higher order BCH can scale rapidly,    so keep this in mind when balancing the need for accuracy and speed of computation.\n",
    "- `truncation_threshold`: This argument allows you to specify a minimum threshold (in terms of error generator rate) below which rates are truncated to zero. This can improve performance      by allowing one to skip the computation of terms corresponding to very small corrections.\n",
    "Some interesting emergent behavior starts to occur when we begin to look at higher-order BCH corrections."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b88d003-76d2-4464-9a77-3c4ce6e63745",
   "metadata": {},
   "outputs": [],
   "source": [
    "propagated_errorgen_layer_second_order = errorgen_propagator.propagate_errorgens_bch(c, bch_order=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a60a527-051b-47e3-bbdc-41492bdeb2f1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(propagated_errorgen_layer_second_order)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90409574-98b3-4fec-9b05-9c5adc9764af",
   "metadata": {},
   "source": [
    "Aside from the fact that there are now significantly more terms than was found for the first-order BCH approximation, notice that there are also now emergent second (and higher) order contributions due to C and A error generators which arise from the composition of purely H and S error generators. These additional terms arise from the non-commutivity of the elementary error generators, particularly the non-commutivity of H and S elementary error generators. For more on this phenomenon see [insert paper reference here]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e818c20-a1a7-4b1c-a068-e0ba548614f8",
   "metadata": {},
   "source": [
    "## Approximate Probabilities and Expectation Values\n",
    "Now you have an efficient representation for an approximation to the effective end-of-circuit error generator for your circuit, what can you do with it? In this section we show how to use this sparse representation to efficiently compute corrections to the outcome probability distributions and pauli observable expectation values of noisy clifford circuits."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "405b95e6-15fc-4a2e-84d5-766a0573ee00",
   "metadata": {},
   "source": [
    "We'll start off by demonstrating how to perform strong simulation using the results of error generator propagation to estimate the output probabilities for a desired computational basis state. \n",
    "\n",
    "To do so we'll be making use of the function `approximate_stabilizer_probability` from the `errgenproptools` module. This function takes as input the following arguments:\n",
    "\n",
    "- errorgen_dict : A dictionary of elementary error generator coefficients and their corresponding rates (as outputted, for example, by `propagate_errorgens_bch`.\n",
    "- circuit : The circuit to compute the output probability for. This can by a pyGSTi `Circuit` object, or alternatively a `stim.Tableau`.\n",
    "- desired_bitstring : A string corresponding to the desired computational basis state.\n",
    "- order : Order of the taylor series approximation for the exponentiated error generator to use in computing the approximate output probability. In principle this function can compute       arbitary-order approximation (but practically the cost of the computation scales in the order).\n",
    "- truncation_threshold : As described above, this is a minimum value below which contributions are truncated to zero which can sometimes improve performance by reducing the number of terms   computed with very small overall corrections to the calculated probability. \n",
    "\n",
    "Let's use the results of the application of the second-order BCH approximation above and compute the approximate probability of reading out the all-zeros state from our circuit. For the ideal circuit, the probability of observing the all-zeros state is 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4cc9e0a0-2da5-40fc-8556-4d3272f1b1be",
   "metadata": {},
   "outputs": [],
   "source": [
    "first_order_approximate_prob = eprop.approximate_stabilizer_probability(propagated_errorgen_layer_second_order, c, '0000', order=1)\n",
    "print(first_order_approximate_prob)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce153eee-90c3-4d8e-b25e-8ec849cde6c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "second_order_approximate_prob = eprop.approximate_stabilizer_probability(propagated_errorgen_layer_second_order, c, '0000', order=2)\n",
    "print(second_order_approximate_prob)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "091dea97-7421-45f3-a1ab-0d17f90d1e88",
   "metadata": {},
   "source": [
    "In this few qubit test case we also have the luxury compare this to the results of the (effectively) exact forward simulation for the error model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2978aee0-7447-452e-85b4-41bbbb79f738",
   "metadata": {},
   "outputs": [],
   "source": [
    "exact_probability = error_model.sim.probs(c)['0000']\n",
    "print(exact_probability)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5fb5fb21-bc96-446b-859c-2252c385c55a",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'Absolute Error Approx to Exact (First-order Taylor, Second-order BCH): {abs(exact_probability-first_order_approximate_prob)}')\n",
    "print(f'Absolute Error Approx to Exact (Second-order Taylor, Second-order BCH): {abs(exact_probability-second_order_approximate_prob)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87e43f33-5b0a-436a-8048-7836879fd205",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'Relative Error Approx to Exact (First-order taylor, Second-order BCH): {100*abs(exact_probability-first_order_approximate_prob)/exact_probability}%')\n",
    "print(f'Relative Error Approx to Exact (Second-order taylor, Second-order BCH): {100*abs(exact_probability-second_order_approximate_prob)/exact_probability}%')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6bc34c95-9553-4c82-a80c-539c801b9029",
   "metadata": {},
   "source": [
    "Here we can see that with the combination of second-order BCH and second-order taylor approximations our estimated probability is accurate to well below a 1 percent relative error. By going out to higher-order in either approximation one can achieve even higher levels of accuracy."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f858f5fd-2bb0-4e38-bab6-9b061aa5a273",
   "metadata": {},
   "source": [
    "In addition to strong simulation of the output probabilities of computational basis states, it is also possible to compute approximate values for the expectation values of pauli observables. The main function for doing so is `approximate_stabilizer_pauli_expectation` from the `errgenproptools` module, the signature of which is nearly identical to that of `approximate_stabilizer_probability` described above, except taking instead a desired pauli observable to estimate the expectation value for. Here we'll again use the results of the second-order BCH approximation produced above and look are various order of the taylor series approximation for the pauli expectation value of 'XYZI' (the value for the ideal noise-free circuit is 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa7c8e66-de23-41b5-8ec8-401d6e3c1bf7",
   "metadata": {},
   "outputs": [],
   "source": [
    "first_order_approximate_pauli_expectation = eprop.approximate_stabilizer_pauli_expectation(propagated_errorgen_layer_second_order, c, 'XYZI', order=1)\n",
    "print(first_order_approximate_pauli_expectation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b6bd6fb-f7e6-482a-97c4-20bcba5749a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "second_order_approximate_pauli_expectation = eprop.approximate_stabilizer_pauli_expectation(propagated_errorgen_layer_second_order, c, 'XYZI', order=2)\n",
    "print(second_order_approximate_pauli_expectation)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c1efde8-7cf4-430a-94a6-376bdf991e67",
   "metadata": {},
   "source": [
    "There aren't existing built-in functions in pyGSTi for outputing exact pauli expectation values handy, but we can write a short helper function for computing these for the sake of comparison with our above results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17731ccb-2cd7-4d7e-8ed8-f21b199d95a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pygsti.tools.basistools import change_basis\n",
    "import numpy as np\n",
    "from pygsti.baseobjs import Label\n",
    "def pauli_expectation_exact(error_propagator, target_model, circuit, pauli):\n",
    "    #get the eoc error channel, and the process matrix for the ideal circuit:\n",
    "    eoc_channel = error_propagator.eoc_error_channel(circuit, include_spam=True)\n",
    "    ideal_channel = target_model.sim.product(circuit)\n",
    "    #also get the ideal state prep and povm:\n",
    "    ideal_prep = target_model.circuit_layer_operator(Label('rho0'), typ='prep').copy()\n",
    "    \n",
    "    #finally need the superoperator for the selected pauli.\n",
    "    pauli_unitary = pauli.to_unitary_matrix(endian='big')\n",
    "    #flatten this row-wise\n",
    "    pauli_vec = np.ravel(pauli_unitary)\n",
    "    pauli_vec.reshape((len(pauli_vec),1))\n",
    "    #put this in pp basis (since these are paulis themselves I could just read this off directly).\n",
    "    pauli_vec = change_basis(pauli_vec, 'std', 'pp')\n",
    "    #print(pauli_vec)\n",
    "    dense_prep = ideal_prep.to_dense().copy()\n",
    "    expectation = np.linalg.multi_dot([pauli_vec.reshape((1,len(pauli_vec))), eoc_channel, ideal_channel, dense_prep.reshape((len(dense_prep),1))]).item()\n",
    "    return expectation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ca1ba09-e416-475d-bb57-58eb28db08a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "exact_pauli_expectation = pauli_expectation_exact(errorgen_propagator, target_model, c, stim.PauliString('XYZI'))\n",
    "print(exact_pauli_expectation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "407e0d72-ced0-4105-8a4f-463085406f18",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'Absolute Error Approx to Exact (First-order Taylor, Second-order BCH): {abs(exact_pauli_expectation-first_order_approximate_pauli_expectation)}')\n",
    "print(f'Absolute Error Approx to Exact (Second-order Taylor, Second-order BCH): {abs(exact_pauli_expectation-second_order_approximate_pauli_expectation)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf63b01d-244c-4943-b477-576b1be496f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'Relative Error Approx to Exact (First-order taylor, Second-order BCH): {100*abs(exact_pauli_expectation-first_order_approximate_pauli_expectation)/exact_pauli_expectation}%')\n",
    "print(f'Relative Error Approx to Exact (Second-order taylor, Second-order BCH): {100*abs(exact_pauli_expectation-second_order_approximate_pauli_expectation)/exact_pauli_expectation}%')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da73f141-80a7-4383-94ae-a530273f3e3d",
   "metadata": {},
   "source": [
    "In this case even with the first-order taylor approximation together with the second-order BCH approximation the relative error to the exact expecation value is roughly half a percent, dropping to below a tenth of a percent when we go up to the second order taylor approximation. As before, by going out to higher-order in either approximation one can achieve even higher levels of accuracy."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a4f38cd-4625-4ecf-9f69-3046d4a45ebd",
   "metadata": {},
   "source": [
    "## Other Helpful Utilities:\n",
    "In this section we'll highlight a few additional utilities within the error generator propagation related modules which are often useful (some of these you may have even seen us use above!).\n",
    "\n",
    "We'll specifically cover:\n",
    "- `eoc_error_channel`\n",
    "- `errorgen_layer_dict_to_errorgen`\n",
    "- `approximate_stabilizer_probabilities`\n",
    "- `error_generator_commutator`\n",
    "- `error_generator_composition`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8ab436d-9b2c-4125-9fb8-69e281225308",
   "metadata": {},
   "source": [
    "#### `eoc_error_channel` : \n",
    "This method provides a simple single function call for generating a dense representation of the end-of-circuit error channel (i.e. the exponentiated end-of-circuit error generator). This can be useful in few-qubit testing, but obviously doesn't not scale beyond a few qubits. This end-of-circuit error channel can be produced either exactly or without the BCH approximation. In the former case this is acheived by exponentiating and multiplying together all of the propagated error generator layers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed42df89-31b9-48af-9f3d-1baf538dc64b",
   "metadata": {},
   "outputs": [],
   "source": [
    "dense_end_of_circuit_channel_exact = errorgen_propagator.eoc_error_channel(c, use_bch=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73cd2961-5280-4d38-8f55-49b6fc4fb7b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "dense_end_of_circuit_channel_first_order_BCH = errorgen_propagator.eoc_error_channel(c, use_bch=True, bch_kwargs={'bch_order':1})\n",
    "dense_end_of_circuit_channel_second_order_BCH = errorgen_propagator.eoc_error_channel(c, use_bch=True, bch_kwargs={'bch_order':2})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3a2d84e-2266-4337-9be7-08f7047141b4",
   "metadata": {},
   "source": [
    "This can be useful in testing settings, for example, where we can use these as yet another way to measure the accuracy of our approximation methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38926ad9-09e1-4174-b845-47bc2aad305e",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'Frobenius norm between exact and 1st-order BCH EOC channels: {np.linalg.norm(dense_end_of_circuit_channel_exact-dense_end_of_circuit_channel_first_order_BCH)}')\n",
    "print(f'Frobenius norm between exact and 2nd-order BCH EOC channels: {np.linalg.norm(dense_end_of_circuit_channel_exact-dense_end_of_circuit_channel_second_order_BCH)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1649c1e7-29a3-492d-b4fb-0cfe2b491c7e",
   "metadata": {},
   "source": [
    "#### `errorgen_layer_dict_to_errorgen`\n",
    "Throughout the error generator propagation framework we generate a lot of sparse error generator representations in terms of dictionaries of elementary error generator coefficients and corresponding rates. For testing purposes (with just a few qubits, this obviously does not scale) it is often useful to convert these into a dense representation as a numpy array. This method helps do so in just a single line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51361d27-6f0f-4765-8b35-81d1f2f62362",
   "metadata": {},
   "outputs": [],
   "source": [
    "dense_end_of_circuit_errorgen_first_order_BCH = errorgen_propagator.errorgen_layer_dict_to_errorgen(propagated_errorgen_layer_first_order)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cedfee7-22db-42a5-98e1-950fdef54149",
   "metadata": {},
   "source": [
    "By default this returns the error generator in the normalized pauli-product basis, but this can be changed using the optional `mx_basis` kwarg.\n",
    "Note: There is another function called `errorgen_layer_to_matrix` available in the `errgenproptools` module with similar functionality to this method, but with a somewhat different interface. That function can be particularly useful in situations where you may want to compute a lot of dense error generator representations from the outputs of the error generator propagation framework, so check out the documentation of that function for more."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2f914cc-2c69-4e80-ab99-9be1fe0b3a9a",
   "metadata": {},
   "source": [
    "#### `approximate_stabilizer_probabilities`\n",
    "This one is straightforward. Above we showed the use of the function `approximate_stabilizer_probability` from the `errgenproptools` module for calculating approximate output probabilities for a given computational bitstring. If you happen to want *all* of the bit string probabilities you can save yourself a for loop by using the function `approximate_stabilizer_probabilities` from this module instead!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21fcd7cc-6799-4aca-b524-69868d0a0169",
   "metadata": {},
   "outputs": [],
   "source": [
    "approximate_probabilities = eprop.approximate_stabilizer_probabilities(propagated_errorgen_layer_first_order, c, order=1)\n",
    "print(approximate_probabilities)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "464e2b37-3d00-4503-9a4d-40f60aa4c3a3",
   "metadata": {},
   "source": [
    "Note the returned values are given in right-LSB convention (i.e. '0000' -> '0001' ->'0010', etc.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3984dac3-6114-4f4e-80e0-4ba31a79886a",
   "metadata": {},
   "source": [
    "#### `error_generator_commutator` and `error_generator_composition`\n",
    "These two functions from the `errgenproptools` module return the result of analytically computing the commutator and composition of two elementary error generators, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83ba69f1-bbd0-4133-8ef3-abd9dc020888",
   "metadata": {},
   "outputs": [],
   "source": [
    "errorgen_1 = _LSE('H', [stim.PauliString('X')])\n",
    "errorgen_2 = _LSE('S', [stim.PauliString('Z')])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ecd3b00-4a7b-4c8e-8a5f-bdba13565fb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(eprop.error_generator_commutator(errorgen_1, errorgen_2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9862ff3-8717-48e7-b0a8-ec3a2d07c974",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(eprop.error_generator_composition(errorgen_1, errorgen_2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a020355-1796-4ddd-ad85-178b37266e35",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(eprop.error_generator_composition(errorgen_1, errorgen_1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9723d08d-09ef-4ea7-ad0e-5624fa6f8501",
   "metadata": {},
   "source": [
    "Both of these methods return their output as a list of two-element tuples. This list is a specification for the linear combination of elementary error generator coefficients corresponding to the commutator or composition of the two input elementary error generators. (First tuple element is an elementary error generator in the linear combination, and the second element is the coefficient of that elementary error generator in the linear combination).\n",
    "\n",
    "In the examples above we can see that the commutator of the specified H and S error generators gives rise to a pauli-correlation (C) error generator. This could potentially give rise to emergent C error generators when applying second-or-higher order BCH approximations for the effective end-of-circuit error generator, for example. Likewise the composition of these to error generators is a linear combination of a C error generator and an H error generator. And finally we see that squaring an H error generator (composing it with itself) gives rise to a pauli-stochastic (S) error generator."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdfce65a-f619-4ef8-b52d-dff5d02a314f",
   "metadata": {},
   "source": [
    "There's a whole bunch of other functionality and utilities available, particularly in the `errgenproptools` module which have not been covered in this tutorial, so please check out the documentation for additional capabilities!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f305946-e23d-4572-909e-6389dfa6a26b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}