From b1dae6647fff8fcda725df6c7cd2340c50b5b42e Mon Sep 17 00:00:00 2001 From: Pablo Andres-Martinez <104848389+PabloAndresCQ@users.noreply.github.com> Date: Mon, 23 Sep 2024 16:15:44 +0100 Subject: [PATCH] Feature/classical conditionals (#152) # Description Adding support so that `simulate` can: - run circuits with classical logic (#151) - run circuits with a single qubit (#156) Also included some bugfixes of MPSxMPO (I had missed calling `_flush()` in the newer functions). # Related issues #151 #156 # Checklist - [x] I have run the tests on a device with GPUs. - [x] I have performed a self-review of my code. - [x] I have commented hard-to-understand parts of my code. - [x] I have made corresponding changes to the public API documentation. - [x] I have added tests that prove my fix is effective or that my feature works. - [x] I have updated the changelog with any user-facing changes. --- docs/changelog.rst | 14 + docs/modules/structured_state.rst | 2 +- examples/mps_tutorial.ipynb | 2 +- examples/python/mps_tutorial.py | 39 +- .../cutensornet/structured_state/classical.py | 135 +++++ .../cutensornet/structured_state/general.py | 99 +++- .../cutensornet/structured_state/mps.py | 65 +-- .../cutensornet/structured_state/mps_mpo.py | 58 +- .../structured_state/simulation.py | 73 ++- .../cutensornet/structured_state/ttn.py | 45 +- tests/conftest.py | 13 + tests/test_structured_state.py | 11 + tests/test_structured_state_conditionals.py | 505 ++++++++++++++++++ 13 files changed, 935 insertions(+), 126 deletions(-) create mode 100644 pytket/extensions/cutensornet/structured_state/classical.py create mode 100644 tests/test_structured_state_conditionals.py diff --git a/docs/changelog.rst b/docs/changelog.rst index c9a51d94..61c44d28 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -1,6 +1,20 @@ Changelog ~~~~~~~~~ +Unreleased +---------- + +* API breaking changes + * Removed ``use_kahypar`` option from ``Config``. It can still be set via the ``simulate`` option ``compilation_params``. + +* New feature: ``simulate`` now accepts pytket circuits with ``Measure``, ``Reset``, ``Conditional``, ``ClassicalExpBox`` and more classical operations. You can now retrieve classical bit values using ``get_bits``. +* When calling ``simulate``, the gates on the circuit are no longer sorted by default. Use ``compilation_params["sort_gates"] = True`` to recover this behaviour, which is now deprecated. +* ``StructuredState`` now supports simulation of single qubit circuits. +* Some bugfixes on MPSxMPO relating to measurement and relabelling qubits. The bug was caused due to these functions not guaranteeing the MPO was applied before their action. +* Documentation fixes: + * ``apply_qubit_relabelling`` now appears in the documentation. + * ``add_qubit`` removed from documentation of MPSxMPO, since it is not supported. + 0.7.1 (July 2024) ----------------- diff --git a/docs/modules/structured_state.rst b/docs/modules/structured_state.rst index 38d22104..8fd0506a 100644 --- a/docs/modules/structured_state.rst +++ b/docs/modules/structured_state.rst @@ -26,6 +26,7 @@ Classes .. automethod:: apply_gate .. automethod:: apply_unitary .. automethod:: apply_scalar + .. automethod:: apply_qubit_relabelling .. automethod:: vdot .. automethod:: sample .. automethod:: measure @@ -52,7 +53,6 @@ Classes .. autoclass:: pytket.extensions.cutensornet.structured_state.MPSxMPO() .. automethod:: __init__ - .. automethod:: add_qubit Miscellaneous diff --git a/examples/mps_tutorial.ipynb b/examples/mps_tutorial.ipynb index db9aacf8..1735b1f9 100644 --- a/examples/mps_tutorial.ipynb +++ b/examples/mps_tutorial.ipynb @@ -1 +1 @@ -{"cells": [{"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["import numpy as np\n", "from time import time\n", "import matplotlib.pyplot as plt\n", "from pytket import Circuit\n", "from pytket.circuit.display import render_circuit_jupyter"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["from pytket.extensions.cutensornet.structured_state import (\n", " CuTensorNetHandle,\n", " Config,\n", " SimulationAlgorithm,\n", " simulate,\n", " prepare_circuit_mps,\n", ")"]}, {"cell_type": "markdown", "metadata": {}, "source": ["# Introduction
\n", "This notebook provides examples of the usage of the MPS functionalities of `pytket_cutensornet`. For more information, see the docs at https://tket.quantinuum.com/extensions/pytket-cutensornet/api/index.html.
\n", "A Matrix Product State (MPS) represents a state on `n` qubits as a list of `n` tensors connected in a line as show below:
\n", "![MPS](images/mps.png)
\n", "Each of these circles corresponds to a tensor. We refer to each leg of a tensor as a *bond* and the number of bonds a tensor has is its *rank*. In code, a tensor is just a multidimensional array:
\n", "
\n", "```tensor[i][j][k] = v```
\n", "
\n", "In the case above, we are assigning an entry value `v` of a rank-3 tensor (one `[ ]` coordinate per bond). Each bond allows a different number of values for its indices; for instance `0 <= i < 4` would mean that the first bond of our tensor can take up to four different indices; we refer to this as the *dimension* of the bond. We refer to the bonds connecting different tensors in the MPS as *virtual bonds*; the maximum allowed value for the dimension of virtual bonds is often denoted by the greek letter `chi`. The open bonds are known as *physical bonds* and, in our case, each will correspond to a qubit; hence, they have dimension `2` -- the dimension of the vector space of a single qubit.
\n", "In essence, whenever we want to apply a gate to certain qubit we will connect a tensor (matrix) representing the gate to the corresponding physical bond and *contract* the network back to an MPS form (tensor contraction is a generalisation of matrix multiplication to multidimensional arrays). Whenever a two-qubit gate is applied, the entanglement information after contraction will be kept in the degrees of freedom of the virtual bonds. As such, the dimension of the virtual bonds will generally increase exponentially as we apply entangling gates, leading to large memory footprints of the tensors and, consequently, long runtime for tensor contraction. We provide functionalities to limit the growth of the dimension of the virtual bonds, keeping resource consumption in check. Read the *Approximate simulation* section on this notebook to learn more.
\n", "**References**: To read more about MPS we recommend the following papers.
\n", "* For an introduction to MPS and its canonical form: https://arxiv.org/abs/1901.05824.
\n", "* For a description of the `MPSxGate` algorithm we provide: https://arxiv.org/abs/2002.07730.
\n", "* For a description of the `MPSxMPO` algorithm we provide: https://arxiv.org/abs/2207.05612.
\n", "* For insights on the reationship between truncation error and the error model in a quantum computer: https://arxiv.org/abs/2004.02388"]}, {"cell_type": "markdown", "metadata": {}, "source": ["# Basic functionality and exact simulation
\n", "Here we show an example of the basic use of our MPS methods. We first generate a simple `pytket` circuit to be simulated."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["my_circ = Circuit(5)\n", "my_circ.CX(3, 4)\n", "my_circ.H(2)\n", "my_circ.CZ(0, 1)\n", "my_circ.ZZPhase(0.1, 4, 3)\n", "my_circ.TK2(0.3, 0.5, 0.7, 2, 1)\n", "my_circ.Ry(0.2, 0)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["render_circuit_jupyter(my_circ)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["For **exact** simulation, simply call the `simulate` function on the circuit and choose a contraction algorithm. To learn more about the contraction algorithms we provide see the *Contraction algorithms* section of this notebook. You will also need to provide a configuration, the default one is provided by `Config()`. Custom settings of `Config` are discussed in the *Approximate simulation* section.
\n", "**NOTE**: whenever you wish to generate an `MPS` object or execute calculations on it you must do so within a `with CuTensorNetHandle() as libhandle:` block; this will initialise the cuTensorNetwork library for you, and destroy its handles at the end of the `with` block. You will need to pass the `libhandle` to the `MPS` object via the method that generates it (in the snippet below, `simulate`), or if already initialised, pass it via the `update_libhandle` method.
\n", "Due to the nature of Jupyter notebooks, we will be starting most of these cells with a `with CuTensorNetHandle() as libhandle:`. However, in a standard script, all of these cells would be grouped together and a single `with CuTensorNetHandle() as libhandle:` statement would be necessary at the beginning of the script."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["with CuTensorNetHandle() as libhandle:\n", " my_mps = simulate(libhandle, my_circ, SimulationAlgorithm.MPSxGate, Config())"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Notice that `my_circ` uses a rich gateset -- in fact, every single-qubit and two-qubit gate supported by `pytket` can be used in our MPS approaches. Gates acting on more than two qubits are not currently supported."]}, {"cell_type": "markdown", "metadata": {}, "source": ["The output of `simulate` is an `MPS` object encoding the output state of the circuit.
\n", "### Obtain an amplitude from an MPS
\n", "Let's first see how to get the amplitude of the state `|10100>` from the output of the previous circuit."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["state = int(\"10100\", 2)\n", "with CuTensorNetHandle() as libhandle:\n", " my_mps.update_libhandle(libhandle)\n", " amplitude = my_mps.get_amplitude(state)\n", "print(amplitude)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Since this is a very small circuit, we can use `pytket`'s state vector simulator capabilities to verify that the state is correct by checking the amplitude of each of the computational states."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["state_vector = my_circ.get_statevector()\n", "n_qubits = len(my_circ.qubits)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["correct_amplitude = [False] * (2**n_qubits)\n", "with CuTensorNetHandle() as libhandle:\n", " my_mps.update_libhandle(libhandle)\n", " for i in range(2**n_qubits):\n", " correct_amplitude[i] = np.isclose(state_vector[i], my_mps.get_amplitude(i))"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["print(\"Are all amplitudes correct?\")\n", "print(all(correct_amplitude))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Sampling from an MPS
\n", "We can also sample from the output state of a circuit by calling `my_mps.sample`, where `my_mps` is the outcome of simulating the circuit."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["n_samples = 100\n", "n_qubits = len(my_circ.qubits)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Initialise the sample counter"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["sample_count = [0 for _ in range(2**n_qubits)]"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["with CuTensorNetHandle() as libhandle:\n", " my_mps.update_libhandle(libhandle)\n", " for _ in range(n_samples):\n", " # Draw a sample\n", " qubit_outcomes = my_mps.sample()\n", " # Convert qubit outcomes to bitstring\n", " bitstring = \"\".join(str(qubit_outcomes[q]) for q in my_circ.qubits)\n", " # Convert bitstring to int\n", " outcome = int(bitstring, 2)\n", " # Update the sample dictionary\n", " sample_count[outcome] += 1"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Calculate the theoretical number of samples per bitstring"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["expected_count = [n_samples * abs(state_vector[i]) ** 2 for i in range(2**n_qubits)]"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Plot a comparison of theory vs sampled"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["plt.scatter(range(2**n_qubits), expected_count, label=\"Theory\")\n", "plt.scatter(range(2**n_qubits), sample_count, label=\"Experiment\", marker=\"x\")\n", "plt.xlabel(\"Basis states\")\n", "plt.ylabel(\"Samples\")\n", "plt.legend()\n", "plt.show()"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We also provide methods to apply mid-circuit measurements via `my_mps.measure(qubits)` and postselection via `my_mps.postselect(qubit_outcomes)`. Their use is similar to that of `my_mps.sample()` shown above.
\n", "**Note:** whereas `my_mps.sample()` does *not* change the state of the MPS, `my_mps.measure(qubits)` and `my_mps.postselect(qubit_outcomes)` do change it, projecting the state to the resulting outcome and removing the measured qubits."]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Inner products
\n", "Using `vdot` you can obtain the inner product of two states in MPS form. This method does not change the internal data of neither of the MPS. Moreover, it can be used on the same `MPS` object for both inputs, yielding the squared norm of the state."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["with CuTensorNetHandle() as libhandle:\n", " my_mps.update_libhandle(libhandle)\n", " norm_sq = my_mps.vdot(my_mps)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["print(\"As expected, the squared norm of a state is 1\")\n", "print(np.isclose(norm_sq, 1))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Let's come up with another circuit on the same qubits and apply an inner product between the two `MPS` objects."]}, {"cell_type": "markdown", "metadata": {}, "source": ["Generate circuits"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["other_circ = Circuit(5)\n", "other_circ.H(3)\n", "other_circ.CZ(3, 4)\n", "other_circ.XXPhase(0.3, 1, 2)\n", "other_circ.Ry(0.7, 3)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Simulate them"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["with CuTensorNetHandle() as libhandle:\n", " other_mps = simulate(libhandle, other_circ, SimulationAlgorithm.MPSxGate, Config())"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Let's calculate the inner product and check that it agrees with `pytket`'s state vector based computation."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["with CuTensorNetHandle() as libhandle:\n", " my_mps.update_libhandle(libhandle)\n", " inner_product = my_mps.vdot(other_mps)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["my_state = my_circ.get_statevector()\n", "other_state = other_circ.get_statevector()"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["print(\"Is the inner product correct?\")\n", "print(np.isclose(np.vdot(my_state, other_state), inner_product))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Two-qubit gates acting on non-adjacent qubits
\n", "Standard MPS algorithms only support simulation of two-qubit gates acting on neighbour qubits. In our implementation, however, two-qubit gates between arbitrary qubits may be applied, as shown below."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["circ = Circuit(5)\n", "circ.H(1)\n", "circ.ZZPhase(0.3, 1, 3)\n", "circ.CX(0, 2)\n", "circ.Ry(0.8, 4)\n", "circ.CZ(3, 4)\n", "circ.XXPhase(0.7, 1, 2)\n", "circ.TK2(0.1, 0.2, 0.4, 1, 4)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["render_circuit_jupyter(circ)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["with CuTensorNetHandle() as libhandle:\n", " mps = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, Config())\n", " print(\"Did simulation succeed?\")\n", " print(mps.is_valid())"]}, {"cell_type": "markdown", "metadata": {}, "source": ["**Note**: Even though two-qubit gates on non-adjacent qubits are simulable, the overhead on these is considerably larger than simulating gates on adjacent qubits. As a rule of thumb if the two qubits are `n` positions apart, the overhead is upper bounded by the cost of simulating `n-1` additional SWAP gates to move the leftmost qubit near the rightmost. In reality, the implementation we use is more nuanced than just applying SWAP gates, and the qubits don't actually change position.
\n", "When circuits are shallow, using our approach to simulate long-distance two-qubit gates is advantageous. In the case of deep circuits with many long-distance gates, it is sometimes beneficial to use TKET routing on the circuit, explicitly adding SWAP gates so that all two-qubit gates act on nearest neighbour qubits. Users may do this by calling `prepare_circuit_mps`, which is a wrapper of the corresponding TKET routing pass."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["prepared_circ, qubit_map = prepare_circuit_mps(circ)\n", "render_circuit_jupyter(prepared_circ)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["The circuit can now be simulated as usual."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["with CuTensorNetHandle() as libhandle:\n", " mps = simulate(libhandle, prepared_circ, SimulationAlgorithm.MPSxGate, Config())\n", " print(\"Did simulation succeed?\")\n", " print(mps.is_valid())"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Notice that the qubits in the `prepared_circ` were renamed when applying `prepare_circuit_mps`. Implicit SWAPs may have been added to the circuit, meaning that the logical qubit held at the ``node[i]`` qubit at the beginning of the circuit may differ from the one it holds at the end; this information is captured by the `qubit_map` output. We recommend applying ``apply_qubit_relabelling`` on the MPS after simulation, relabelling the qubits according to these implicit SWAPs."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["print(qubit_map)\n", "mps.apply_qubit_relabelling(qubit_map)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["# Approximate simulation
\n", "We provide two policies for approximate simulation; these are supported by both of our current MPS contraction algorithms:
\n", "* Bound the maximum value of the virtual bond dimension `chi`. If a bond dimension would increase past that point, we *truncate* (i.e. discard) the degrees of freedom that contribute the least to the state description. We can keep track of a lower bound of the error that this truncation causes.
\n", "* Provide a value for acceptable two-qubit gate fidelity `truncation_fidelity`. After each two-qubit gate we truncate the dimension of virtual bonds as much as we can while guaranteeing the target gate fidelity. The more fidelity you require, the longer it will take to simulate. **Note**: this is *not* the final fidelity of the output state, but the fidelity per gate.
\n", "Values for `chi` and `truncation_fidelity` can be set via `Config`. To showcase approximate simulation, let's define a circuit where exact MPS contraction starts struggling."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["def random_line_circuit(n_qubits: int, layers: int) -> Circuit:\n", " \"\"\"Random circuit with line connectivity.\"\"\"\n", " c = Circuit(n_qubits)\n", " for i in range(layers):\n", " # Layer of TK1 gates\n", " for q in range(n_qubits):\n", " c.TK1(np.random.rand(), np.random.rand(), np.random.rand(), q)\n\n", " # Layer of CX gates\n", " offset = np.mod(i, 2) # Even layers connect (q0,q1), odd (q1,q2)\n", " qubit_pairs = [\n", " [c.qubits[i], c.qubits[i + 1]] for i in range(offset, n_qubits - 1, 2)\n", " ]\n", " # Direction of each CX gate is random\n", " for pair in qubit_pairs:\n", " np.random.shuffle(pair)\n", " for pair in qubit_pairs:\n", " c.CX(pair[0], pair[1])\n", " return c"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["circuit = random_line_circuit(n_qubits=20, layers=20)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["For exact contraction, `chi` must be allowed to be up to `2**(n_qubits // 2)`, meaning that if we set `n_qubits = 20` it would require `chi = 1024`; already too much for this particular circuit to be simulated in a gaming laptop using the current implementation. Instead, let's bound `chi` to a maximum of `16`. Doing so results in faster runtime, at the expense of losing output state fidelity."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["start = time()\n", "with CuTensorNetHandle() as libhandle:\n", " config = Config(chi=16)\n", " bound_chi_mps = simulate(libhandle, circuit, SimulationAlgorithm.MPSxGate, config)\n", "end = time()\n", "print(\"Time taken by approximate contraction with bound chi:\")\n", "print(f\"{round(end-start,2)} seconds\")\n", "print(\"\\nLower bound of the fidelity:\")\n", "print(round(bound_chi_mps.fidelity, 4))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Alternatively, we can fix `truncation_fidelity` and let `chi` increase as necessary to satisfy it."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["start = time()\n", "with CuTensorNetHandle() as libhandle:\n", " config = Config(truncation_fidelity=0.999)\n", " fixed_fidelity_mps = simulate(\n", " libhandle, circuit, SimulationAlgorithm.MPSxGate, config\n", " )\n", "end = time()\n", "print(\"Time taken by approximate contraction with fixed truncation fidelity:\")\n", "print(f\"{round(end-start,2)} seconds\")\n", "print(\"\\nLower bound of the fidelity:\")\n", "print(round(fixed_fidelity_mps.fidelity, 4))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["# Contraction algorithms"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We currently offer two MPS-based simulation algorithms:
\n", "* **MPSxGate**: Apply gates one by one to the MPS, canonicalising the MPS and truncating when necessary. In particular, we implemented the algorithm from the following paper: https://arxiv.org/abs/2002.07730.
\n", "* **MPSxMPO**: Maintain two MPS copies of the state as it evolves, one updated eagerly using the **MPSxGate** method and the other updated in batches of up to `k` layers of two-qubit gates. Whenever the second MPS is updated, both copies are synchronised and an optimisation algorithm is applied to increase the fidelity of the state. This algorithm is often referred to as DMRG-like simulation. In particular, we implemented the algorithm from the following paper: https://arxiv.org/abs/2207.05612.
\n", "The `MPSxGate` algorithm is the one we have been using for all of the examples above. In comparison, the `MPSxMPO` algorithm provides the user with two new parameters to tune:
\n", "* **k**: The maximum number of layers the MPO is allowed to have before being contracted. Increasing this might increase fidelity, but it will also increase resource requirements exponentially. Default value is `4`.
\n", "* **optim_delta**: Stopping criteria for the optimisation when contracting the `k` layers of MPO. Stops when the increase of fidelity between iterations is smaller than `optim_delta`. Default value is `1e-5`.
\n", "Both `k` and `optim_delta` can be set via `Config`. Below we compare `MPSxGate` versus `MPSxMPO` with default parameters and `MPSxMPO` with more resource-hungry parameters. The circuit used is the same as in the previous section."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["start = time()\n", "with CuTensorNetHandle() as libhandle:\n", " config = Config(chi=16)\n", " fixed_fidelity_mps = simulate(\n", " libhandle, circuit, SimulationAlgorithm.MPSxGate, config\n", " )\n", "end = time()\n", "print(\"MPSxGate\")\n", "print(f\"\\tTime taken: {round(end-start,2)} seconds\")\n", "print(f\"\\tLower bound of the fidelity: {round(fixed_fidelity_mps.fidelity, 4)}\")"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["start = time()\n", "with CuTensorNetHandle() as libhandle:\n", " config = Config(chi=16)\n", " fixed_fidelity_mps = simulate(\n", " libhandle, circuit, SimulationAlgorithm.MPSxMPO, config\n", " )\n", "end = time()\n", "print(\"MPSxMPO, default parameters\")\n", "print(f\"\\tTime taken: {round(end-start,2)} seconds\")\n", "print(f\"\\tLower bound of the fidelity: {round(fixed_fidelity_mps.fidelity, 4)}\")"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["start = time()\n", "with CuTensorNetHandle() as libhandle:\n", " config = Config(k=8, optim_delta=1e-15, chi=16)\n", " fixed_fidelity_mps = simulate(\n", " libhandle, circuit, SimulationAlgorithm.MPSxMPO, config\n", " )\n", "end = time()\n", "print(\"MPSxMPO, custom parameters\")\n", "print(f\"\\tTime taken: {round(end-start,2)} seconds\")\n", "print(f\"\\tLower bound of the fidelity: {round(fixed_fidelity_mps.fidelity, 4)}\")"]}, {"cell_type": "markdown", "metadata": {}, "source": ["**Note**: `MPSxMPO` also admits truncation policy in terms of `truncation_fidelity` instead of `chi`."]}, {"cell_type": "markdown", "metadata": {}, "source": ["# Using the logger"]}, {"cell_type": "markdown", "metadata": {}, "source": ["You can request a verbose log to be produced during simulation, by assigning the `loglevel` argument when creating a `Config` instance. Currently, two log levels are supported (other than default, which is silent):
\n", "- `logging.INFO` will print information about progress percent, memory currently occupied by the MPS and current fidelity. Additionally, some high level information of the current stage of the simulation is provided, such as when `MPSxMPO` is applying optimisation sweeps.
\n", "- `logging.DEBUG` provides all of the messages from the loglevel above plus detailed information of the current operation being carried out and the values of important variables.
\n", "**Note**: Due to technical issues with the `logging` module and Jupyter notebooks we need to reload the `logging` module. When working with python scripts and command line, just doing `import logging` is enough."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["from importlib import reload # Not needed in Python 2\n", "import logging"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["reload(logging)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["An example of the use of `logging.INFO` is provided below."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["with CuTensorNetHandle() as libhandle:\n", " config = Config(truncation_fidelity=0.999, loglevel=logging.INFO)\n", " simulate(libhandle, circuit, SimulationAlgorithm.MPSxMPO, config)"]}], "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} \ No newline at end of file +{"cells": [{"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["import numpy as np\n", "from time import time\n", "import matplotlib.pyplot as plt\n", "from pytket import Circuit, OpType\n", "from pytket.circuit.display import render_circuit_jupyter"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["from pytket.extensions.cutensornet.structured_state import (\n", " CuTensorNetHandle,\n", " Config,\n", " SimulationAlgorithm,\n", " simulate,\n", " prepare_circuit_mps,\n", ")"]}, {"cell_type": "markdown", "metadata": {}, "source": ["# Introduction
\n", "This notebook provides examples of the usage of the MPS functionalities of `pytket_cutensornet`. For more information, see the docs at https://tket.quantinuum.com/extensions/pytket-cutensornet/api/index.html.
\n", "A Matrix Product State (MPS) represents a state on `n` qubits as a list of `n` tensors connected in a line as show below:
\n", "![MPS](images/mps.png)
\n", "Each of these circles corresponds to a tensor. We refer to each leg of a tensor as a *bond* and the number of bonds a tensor has is its *rank*. In code, a tensor is just a multidimensional array:
\n", "
\n", "```tensor[i][j][k] = v```
\n", "
\n", "In the case above, we are assigning an entry value `v` of a rank-3 tensor (one `[ ]` coordinate per bond). Each bond allows a different number of values for its indices; for instance `0 <= i < 4` would mean that the first bond of our tensor can take up to four different indices; we refer to this as the *dimension* of the bond. We refer to the bonds connecting different tensors in the MPS as *virtual bonds*; the maximum allowed value for the dimension of virtual bonds is often denoted by the greek letter `chi`. The open bonds are known as *physical bonds* and, in our case, each will correspond to a qubit; hence, they have dimension `2` -- the dimension of the vector space of a single qubit.
\n", "In essence, whenever we want to apply a gate to certain qubit we will connect a tensor (matrix) representing the gate to the corresponding physical bond and *contract* the network back to an MPS form (tensor contraction is a generalisation of matrix multiplication to multidimensional arrays). Whenever a two-qubit gate is applied, the entanglement information after contraction will be kept in the degrees of freedom of the virtual bonds. As such, the dimension of the virtual bonds will generally increase exponentially as we apply entangling gates, leading to large memory footprints of the tensors and, consequently, long runtime for tensor contraction. We provide functionalities to limit the growth of the dimension of the virtual bonds, keeping resource consumption in check. Read the *Approximate simulation* section on this notebook to learn more.
\n", "**References**: To read more about MPS we recommend the following papers.
\n", "* For an introduction to MPS and its canonical form: https://arxiv.org/abs/1901.05824.
\n", "* For a description of the `MPSxGate` algorithm we provide: https://arxiv.org/abs/2002.07730.
\n", "* For a description of the `MPSxMPO` algorithm we provide: https://arxiv.org/abs/2207.05612.
\n", "* For insights on the reationship between truncation error and the error model in a quantum computer: https://arxiv.org/abs/2004.02388"]}, {"cell_type": "markdown", "metadata": {}, "source": ["# Basic functionality and exact simulation
\n", "Here we show an example of the basic use of our MPS methods. We first generate a simple `pytket` circuit to be simulated."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["my_circ = Circuit(5)\n", "my_circ.CX(3, 4)\n", "my_circ.H(2)\n", "my_circ.CZ(0, 1)\n", "my_circ.ZZPhase(0.1, 4, 3)\n", "my_circ.TK2(0.3, 0.5, 0.7, 2, 1)\n", "my_circ.Ry(0.2, 0)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["render_circuit_jupyter(my_circ)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["For **exact** simulation, simply call the `simulate` function on the circuit and choose a contraction algorithm. To learn more about the contraction algorithms we provide see the *Contraction algorithms* section of this notebook. You will also need to provide a configuration, the default one is provided by `Config()`. Custom settings of `Config` are discussed in the *Approximate simulation* section.
\n", "**NOTE**: whenever you wish to generate an `MPS` object or execute calculations on it you must do so within a `with CuTensorNetHandle() as libhandle:` block; this will initialise the cuTensorNetwork library for you, and destroy its handles at the end of the `with` block. You will need to pass the `libhandle` to the `MPS` object via the method that generates it (in the snippet below, `simulate`), or if already initialised, pass it via the `update_libhandle` method.
\n", "Due to the nature of Jupyter notebooks, we will be starting most of these cells with a `with CuTensorNetHandle() as libhandle:`. However, in a standard script, all of these cells would be grouped together and a single `with CuTensorNetHandle() as libhandle:` statement would be necessary at the beginning of the script."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["with CuTensorNetHandle() as libhandle:\n", " my_mps = simulate(libhandle, my_circ, SimulationAlgorithm.MPSxGate, Config())"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Notice that `my_circ` uses a rich gateset -- in fact, every single-qubit and two-qubit gate supported by `pytket` can be used in our MPS approaches. Gates acting on more than two qubits are not currently supported."]}, {"cell_type": "markdown", "metadata": {}, "source": ["The output of `simulate` is an `MPS` object encoding the output state of the circuit.
\n", "### Obtain an amplitude from an MPS
\n", "Let's first see how to get the amplitude of the state `|10100>` from the output of the previous circuit."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["state = int(\"10100\", 2)\n", "with CuTensorNetHandle() as libhandle:\n", " my_mps.update_libhandle(libhandle)\n", " amplitude = my_mps.get_amplitude(state)\n", "print(amplitude)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Since this is a very small circuit, we can use `pytket`'s state vector simulator capabilities to verify that the state is correct by checking the amplitude of each of the computational states."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["state_vector = my_circ.get_statevector()\n", "n_qubits = len(my_circ.qubits)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["correct_amplitude = [False] * (2**n_qubits)\n", "with CuTensorNetHandle() as libhandle:\n", " my_mps.update_libhandle(libhandle)\n", " for i in range(2**n_qubits):\n", " correct_amplitude[i] = np.isclose(state_vector[i], my_mps.get_amplitude(i))"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["print(\"Are all amplitudes correct?\")\n", "print(all(correct_amplitude))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Sampling from an MPS
\n", "We can also sample from the output state of a circuit by calling `my_mps.sample`, where `my_mps` is the outcome of simulating the circuit."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["n_samples = 100\n", "n_qubits = len(my_circ.qubits)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Initialise the sample counter"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["sample_count = [0 for _ in range(2**n_qubits)]"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["with CuTensorNetHandle() as libhandle:\n", " my_mps.update_libhandle(libhandle)\n", " for _ in range(n_samples):\n", " # Draw a sample\n", " qubit_outcomes = my_mps.sample()\n", " # Convert qubit outcomes to bitstring\n", " bitstring = \"\".join(str(qubit_outcomes[q]) for q in my_circ.qubits)\n", " # Convert bitstring to int\n", " outcome = int(bitstring, 2)\n", " # Update the sample dictionary\n", " sample_count[outcome] += 1"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Calculate the theoretical number of samples per bitstring"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["expected_count = [n_samples * abs(state_vector[i]) ** 2 for i in range(2**n_qubits)]"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Plot a comparison of theory vs sampled"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["plt.scatter(range(2**n_qubits), expected_count, label=\"Theory\")\n", "plt.scatter(range(2**n_qubits), sample_count, label=\"Experiment\", marker=\"x\")\n", "plt.xlabel(\"Basis states\")\n", "plt.ylabel(\"Samples\")\n", "plt.legend()\n", "plt.show()"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We also provide methods to apply mid-circuit measurements via `my_mps.measure(qubits)` and postselection via `my_mps.postselect(qubit_outcomes)`. Their use is similar to that of `my_mps.sample()` shown above.
\n", "**Note:** whereas `my_mps.sample()` does *not* change the state of the MPS, `my_mps.measure(qubits)` and `my_mps.postselect(qubit_outcomes)` do change it, projecting the state to the resulting outcome and removing the measured qubits."]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Inner products
\n", "Using `vdot` you can obtain the inner product of two states in MPS form. This method does not change the internal data of neither of the MPS. Moreover, it can be used on the same `MPS` object for both inputs, yielding the squared norm of the state."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["with CuTensorNetHandle() as libhandle:\n", " my_mps.update_libhandle(libhandle)\n", " norm_sq = my_mps.vdot(my_mps)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["print(\"As expected, the squared norm of a state is 1\")\n", "print(np.isclose(norm_sq, 1))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Let's come up with another circuit on the same qubits and apply an inner product between the two `MPS` objects."]}, {"cell_type": "markdown", "metadata": {}, "source": ["Generate circuits"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["other_circ = Circuit(5)\n", "other_circ.H(3)\n", "other_circ.CZ(3, 4)\n", "other_circ.XXPhase(0.3, 1, 2)\n", "other_circ.Ry(0.7, 3)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Simulate them"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["with CuTensorNetHandle() as libhandle:\n", " other_mps = simulate(libhandle, other_circ, SimulationAlgorithm.MPSxGate, Config())"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Let's calculate the inner product and check that it agrees with `pytket`'s state vector based computation."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["with CuTensorNetHandle() as libhandle:\n", " my_mps.update_libhandle(libhandle)\n", " inner_product = my_mps.vdot(other_mps)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["my_state = my_circ.get_statevector()\n", "other_state = other_circ.get_statevector()"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["print(\"Is the inner product correct?\")\n", "print(np.isclose(np.vdot(my_state, other_state), inner_product))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Mid-circuit measurements and classical control
\n", "Mid-circuit measurements and classical control is supported (only in `MPSxGate` as of v0.8.0). For instance, we can implement the teleportation protocol on a pytket circuit and simulate it:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["circ = Circuit()\n", "alice = circ.add_q_register(\"alice\", 2)\n", "alice_bits = circ.add_c_register(\"alice_bits\", 2)\n", "bob = circ.add_q_register(\"bob\", 1)\n", "# Initialise Alice's first qubit in some arbitrary state\n", "circ.Rx(0.42, alice[0])\n", "orig_state = circ.get_statevector()\n", "# Create a Bell pair shared between Alice and Bob\n", "circ.H(alice[1]).CX(alice[1], bob[0])\n", "# Apply a Bell measurement on Alice's qubits\n", "circ.CX(alice[0], alice[1]).H(alice[0])\n", "circ.Measure(alice[0], alice_bits[0])\n", "circ.Measure(alice[1], alice_bits[1])\n", "# Apply conditional corrections on Bob's qubits\n", "circ.add_gate(OpType.X, [bob[0]], condition_bits=[alice_bits[1]], condition_value=1)\n", "circ.add_gate(OpType.Z, [bob[0]], condition_bits=[alice_bits[0]], condition_value=1)\n", "# Reset Alice's qubits\n", "circ.add_gate(OpType.Reset, [alice[0]])\n", "circ.add_gate(OpType.Reset, [alice[1]])\n", "# Display the circuit\n", "render_circuit_jupyter(circ)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We can now simulate the circuit and check that the qubit has been successfully teleported."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["print(\n", " f\"Initial state:\\n {np.round(orig_state[0],2)}|00>|0> + {np.round(orig_state[4],2)}|10>|0>\"\n", ")\n", "with CuTensorNetHandle() as libhandle:\n", " state = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, Config())\n", " print(\n", " f\"Teleported state:\\n {np.round(state.get_amplitude(0),2)}|00>|0> + {np.round(state.get_amplitude(1),2)}|00>|1>\"\n", " )\n", " print(f\"Measurement outcomes:\\n {state.get_bits()}\")"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Two-qubit gates acting on non-adjacent qubits
\n", "Standard MPS algorithms only support simulation of two-qubit gates acting on neighbour qubits. In our implementation, however, two-qubit gates between arbitrary qubits may be applied, as shown below."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["circ = Circuit(5)\n", "circ.H(1)\n", "circ.ZZPhase(0.3, 1, 3)\n", "circ.CX(0, 2)\n", "circ.Ry(0.8, 4)\n", "circ.CZ(3, 4)\n", "circ.XXPhase(0.7, 1, 2)\n", "circ.TK2(0.1, 0.2, 0.4, 1, 4)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["render_circuit_jupyter(circ)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["with CuTensorNetHandle() as libhandle:\n", " mps = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, Config())\n", " print(\"Did simulation succeed?\")\n", " print(mps.is_valid())"]}, {"cell_type": "markdown", "metadata": {}, "source": ["**Note**: Even though two-qubit gates on non-adjacent qubits are simulable, the overhead on these is considerably larger than simulating gates on adjacent qubits. As a rule of thumb if the two qubits are `n` positions apart, the overhead is upper bounded by the cost of simulating `n-1` additional SWAP gates to move the leftmost qubit near the rightmost. In reality, the implementation we use is more nuanced than just applying SWAP gates, and the qubits don't actually change position.
\n", "When circuits are shallow, using our approach to simulate long-distance two-qubit gates is advantageous. In the case of deep circuits with many long-distance gates, it is sometimes beneficial to use TKET routing on the circuit, explicitly adding SWAP gates so that all two-qubit gates act on nearest neighbour qubits. Users may do this by calling `prepare_circuit_mps`, which is a wrapper of the corresponding TKET routing pass."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["prepared_circ, qubit_map = prepare_circuit_mps(circ)\n", "render_circuit_jupyter(prepared_circ)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["The circuit can now be simulated as usual."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["with CuTensorNetHandle() as libhandle:\n", " mps = simulate(libhandle, prepared_circ, SimulationAlgorithm.MPSxGate, Config())\n", " print(\"Did simulation succeed?\")\n", " print(mps.is_valid())"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Notice that the qubits in the `prepared_circ` were renamed when applying `prepare_circuit_mps`. Implicit SWAPs may have been added to the circuit, meaning that the logical qubit held at the ``node[i]`` qubit at the beginning of the circuit may differ from the one it holds at the end; this information is captured by the `qubit_map` output. We recommend applying ``apply_qubit_relabelling`` on the MPS after simulation, relabelling the qubits according to these implicit SWAPs."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["print(qubit_map)\n", "mps.apply_qubit_relabelling(qubit_map)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["# Approximate simulation
\n", "We provide two policies for approximate simulation; these are supported by both of our current MPS contraction algorithms:
\n", "* Bound the maximum value of the virtual bond dimension `chi`. If a bond dimension would increase past that point, we *truncate* (i.e. discard) the degrees of freedom that contribute the least to the state description. We can keep track of a lower bound of the error that this truncation causes.
\n", "* Provide a value for acceptable two-qubit gate fidelity `truncation_fidelity`. After each two-qubit gate we truncate the dimension of virtual bonds as much as we can while guaranteeing the target gate fidelity. The more fidelity you require, the longer it will take to simulate. **Note**: this is *not* the final fidelity of the output state, but the fidelity per gate.
\n", "Values for `chi` and `truncation_fidelity` can be set via `Config`. To showcase approximate simulation, let's define a circuit where exact MPS contraction starts struggling."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["def random_line_circuit(n_qubits: int, layers: int) -> Circuit:\n", " \"\"\"Random circuit with line connectivity.\"\"\"\n", " c = Circuit(n_qubits)\n", " for i in range(layers):\n", " # Layer of TK1 gates\n", " for q in range(n_qubits):\n", " c.TK1(np.random.rand(), np.random.rand(), np.random.rand(), q)\n\n", " # Layer of CX gates\n", " offset = np.mod(i, 2) # Even layers connect (q0,q1), odd (q1,q2)\n", " qubit_pairs = [\n", " [c.qubits[i], c.qubits[i + 1]] for i in range(offset, n_qubits - 1, 2)\n", " ]\n", " # Direction of each CX gate is random\n", " for pair in qubit_pairs:\n", " np.random.shuffle(pair)\n", " for pair in qubit_pairs:\n", " c.CX(pair[0], pair[1])\n", " return c"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["circuit = random_line_circuit(n_qubits=20, layers=20)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["For exact contraction, `chi` must be allowed to be up to `2**(n_qubits // 2)`, meaning that if we set `n_qubits = 20` it would require `chi = 1024`; already too much for this particular circuit to be simulated in a gaming laptop using the current implementation. Instead, let's bound `chi` to a maximum of `16`. Doing so results in faster runtime, at the expense of losing output state fidelity."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["start = time()\n", "with CuTensorNetHandle() as libhandle:\n", " config = Config(chi=16)\n", " bound_chi_mps = simulate(libhandle, circuit, SimulationAlgorithm.MPSxGate, config)\n", "end = time()\n", "print(\"Time taken by approximate contraction with bound chi:\")\n", "print(f\"{round(end-start,2)} seconds\")\n", "print(\"\\nLower bound of the fidelity:\")\n", "print(round(bound_chi_mps.fidelity, 4))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Alternatively, we can fix `truncation_fidelity` and let `chi` increase as necessary to satisfy it."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["start = time()\n", "with CuTensorNetHandle() as libhandle:\n", " config = Config(truncation_fidelity=0.999)\n", " fixed_fidelity_mps = simulate(\n", " libhandle, circuit, SimulationAlgorithm.MPSxGate, config\n", " )\n", "end = time()\n", "print(\"Time taken by approximate contraction with fixed truncation fidelity:\")\n", "print(f\"{round(end-start,2)} seconds\")\n", "print(\"\\nLower bound of the fidelity:\")\n", "print(round(fixed_fidelity_mps.fidelity, 4))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["# Contraction algorithms"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We currently offer two MPS-based simulation algorithms:
\n", "* **MPSxGate**: Apply gates one by one to the MPS, canonicalising the MPS and truncating when necessary. In particular, we implemented the algorithm from the following paper: https://arxiv.org/abs/2002.07730.
\n", "* **MPSxMPO**: Maintain two MPS copies of the state as it evolves, one updated eagerly using the **MPSxGate** method and the other updated in batches of up to `k` layers of two-qubit gates. Whenever the second MPS is updated, both copies are synchronised and an optimisation algorithm is applied to increase the fidelity of the state. This algorithm is often referred to as DMRG-like simulation. In particular, we implemented the algorithm from the following paper: https://arxiv.org/abs/2207.05612.
\n", "The `MPSxGate` algorithm is the one we have been using for all of the examples above. In comparison, the `MPSxMPO` algorithm provides the user with two new parameters to tune:
\n", "* **k**: The maximum number of layers the MPO is allowed to have before being contracted. Increasing this might increase fidelity, but it will also increase resource requirements exponentially. Default value is `4`.
\n", "* **optim_delta**: Stopping criteria for the optimisation when contracting the `k` layers of MPO. Stops when the increase of fidelity between iterations is smaller than `optim_delta`. Default value is `1e-5`.
\n", "Both `k` and `optim_delta` can be set via `Config`. Below we compare `MPSxGate` versus `MPSxMPO` with default parameters and `MPSxMPO` with more resource-hungry parameters. The circuit used is the same as in the previous section."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["start = time()\n", "with CuTensorNetHandle() as libhandle:\n", " config = Config(chi=16)\n", " fixed_fidelity_mps = simulate(\n", " libhandle, circuit, SimulationAlgorithm.MPSxGate, config\n", " )\n", "end = time()\n", "print(\"MPSxGate\")\n", "print(f\"\\tTime taken: {round(end-start,2)} seconds\")\n", "print(f\"\\tLower bound of the fidelity: {round(fixed_fidelity_mps.fidelity, 4)}\")"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["start = time()\n", "with CuTensorNetHandle() as libhandle:\n", " config = Config(chi=16)\n", " fixed_fidelity_mps = simulate(\n", " libhandle, circuit, SimulationAlgorithm.MPSxMPO, config\n", " )\n", "end = time()\n", "print(\"MPSxMPO, default parameters\")\n", "print(f\"\\tTime taken: {round(end-start,2)} seconds\")\n", "print(f\"\\tLower bound of the fidelity: {round(fixed_fidelity_mps.fidelity, 4)}\")"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["start = time()\n", "with CuTensorNetHandle() as libhandle:\n", " config = Config(k=8, optim_delta=1e-15, chi=16)\n", " fixed_fidelity_mps = simulate(\n", " libhandle, circuit, SimulationAlgorithm.MPSxMPO, config\n", " )\n", "end = time()\n", "print(\"MPSxMPO, custom parameters\")\n", "print(f\"\\tTime taken: {round(end-start,2)} seconds\")\n", "print(f\"\\tLower bound of the fidelity: {round(fixed_fidelity_mps.fidelity, 4)}\")"]}, {"cell_type": "markdown", "metadata": {}, "source": ["**Note**: `MPSxMPO` also admits truncation policy in terms of `truncation_fidelity` instead of `chi`."]}, {"cell_type": "markdown", "metadata": {}, "source": ["# Using the logger"]}, {"cell_type": "markdown", "metadata": {}, "source": ["You can request a verbose log to be produced during simulation, by assigning the `loglevel` argument when creating a `Config` instance. Currently, two log levels are supported (other than default, which is silent):
\n", "- `logging.INFO` will print information about progress percent, memory currently occupied by the MPS and current fidelity. Additionally, some high level information of the current stage of the simulation is provided, such as when `MPSxMPO` is applying optimisation sweeps.
\n", "- `logging.DEBUG` provides all of the messages from the loglevel above plus detailed information of the current operation being carried out and the values of important variables.
\n", "**Note**: Due to technical issues with the `logging` module and Jupyter notebooks we need to reload the `logging` module. When working with python scripts and command line, just doing `import logging` is enough."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["from importlib import reload # Not needed in Python 2\n", "import logging"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["reload(logging)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["An example of the use of `logging.INFO` is provided below."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["with CuTensorNetHandle() as libhandle:\n", " config = Config(truncation_fidelity=0.999, loglevel=logging.INFO)\n", " simulate(libhandle, circuit, SimulationAlgorithm.MPSxMPO, config)"]}], "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} \ No newline at end of file diff --git a/examples/python/mps_tutorial.py b/examples/python/mps_tutorial.py index 99a39fc3..21a1daca 100644 --- a/examples/python/mps_tutorial.py +++ b/examples/python/mps_tutorial.py @@ -1,7 +1,7 @@ import numpy as np from time import time import matplotlib.pyplot as plt -from pytket import Circuit +from pytket import Circuit, OpType from pytket.circuit.display import render_circuit_jupyter from pytket.extensions.cutensornet.structured_state import ( @@ -145,6 +145,43 @@ print("Is the inner product correct?") print(np.isclose(np.vdot(my_state, other_state), inner_product)) +# ### Mid-circuit measurements and classical control +# Mid-circuit measurements and classical control is supported (only in `MPSxGate` as of v0.8.0). For instance, we can implement the teleportation protocol on a pytket circuit and simulate it: + +circ = Circuit() +alice = circ.add_q_register("alice", 2) +alice_bits = circ.add_c_register("alice_bits", 2) +bob = circ.add_q_register("bob", 1) +# Initialise Alice's first qubit in some arbitrary state +circ.Rx(0.42, alice[0]) +orig_state = circ.get_statevector() +# Create a Bell pair shared between Alice and Bob +circ.H(alice[1]).CX(alice[1], bob[0]) +# Apply a Bell measurement on Alice's qubits +circ.CX(alice[0], alice[1]).H(alice[0]) +circ.Measure(alice[0], alice_bits[0]) +circ.Measure(alice[1], alice_bits[1]) +# Apply conditional corrections on Bob's qubits +circ.add_gate(OpType.X, [bob[0]], condition_bits=[alice_bits[1]], condition_value=1) +circ.add_gate(OpType.Z, [bob[0]], condition_bits=[alice_bits[0]], condition_value=1) +# Reset Alice's qubits +circ.add_gate(OpType.Reset, [alice[0]]) +circ.add_gate(OpType.Reset, [alice[1]]) +# Display the circuit +render_circuit_jupyter(circ) + +# We can now simulate the circuit and check that the qubit has been successfully teleported. + +print( + f"Initial state:\n {np.round(orig_state[0],2)}|00>|0> + {np.round(orig_state[4],2)}|10>|0>" +) +with CuTensorNetHandle() as libhandle: + state = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, Config()) + print( + f"Teleported state:\n {np.round(state.get_amplitude(0),2)}|00>|0> + {np.round(state.get_amplitude(1),2)}|00>|1>" + ) + print(f"Measurement outcomes:\n {state.get_bits()}") + # ### Two-qubit gates acting on non-adjacent qubits # Standard MPS algorithms only support simulation of two-qubit gates acting on neighbour qubits. In our implementation, however, two-qubit gates between arbitrary qubits may be applied, as shown below. diff --git a/pytket/extensions/cutensornet/structured_state/classical.py b/pytket/extensions/cutensornet/structured_state/classical.py new file mode 100644 index 00000000..d49ac4d9 --- /dev/null +++ b/pytket/extensions/cutensornet/structured_state/classical.py @@ -0,0 +1,135 @@ +# Copyright 2019-2024 Quantinuum +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +## +# http://www.apache.org/licenses/LICENSE-2.0 +## +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from typing import Union, Any + +from pytket.circuit import ( + Op, + OpType, + Bit, + BitRegister, + SetBitsOp, + CopyBitsOp, + RangePredicateOp, + ClassicalExpBox, + LogicExp, + BitWiseOp, + RegWiseOp, +) + + +ExtendedLogicExp = Union[LogicExp, Bit, BitRegister, int] + + +def apply_classical_command( + op: Op, bits: list[Bit], args: list[Any], bits_dict: dict[Bit, bool] +) -> None: + """Evaluate classical commands and update the `bits_dict` accordingly.""" + if isinstance(op, SetBitsOp): + for b, v in zip(bits, op.values): + bits_dict[b] = v + + elif isinstance(op, CopyBitsOp): + output_bits = bits + input_bits = args[: len(output_bits)] + for i, o in zip(input_bits, output_bits): + assert isinstance(i, Bit) + bits_dict[o] = bits_dict[i] + + elif isinstance(op, RangePredicateOp): + assert len(bits) == 1 + res_bit = bits[0] + input_bits = args[:-1] + # The input_bits encode a "value" int in little-endian + val = from_little_endian([bits_dict[b] for b in input_bits]) # type: ignore + # Check that the value is in the range + bits_dict[res_bit] = val >= op.lower and val <= op.upper + + elif isinstance(op, ClassicalExpBox): + the_exp = op.get_exp() + result = evaluate_logic_exp(the_exp, bits_dict) + + # The result is an int in little-endian encoding. We update the + # output register accordingly. + for b in bits: + bits_dict[b] = (result % 2) == 1 + result = result >> 1 + assert result == 0 # All bits consumed + + elif op.type == OpType.Barrier: + pass + + else: + raise NotImplementedError(f"Commands of type {op.type} are not supported.") + + +def evaluate_logic_exp(exp: ExtendedLogicExp, bits_dict: dict[Bit, bool]) -> int: + """Recursive evaluation of a LogicExp.""" + + if isinstance(exp, int): + return exp + elif isinstance(exp, Bit): + return 1 if bits_dict[exp] else 0 + elif isinstance(exp, BitRegister): + return from_little_endian([bits_dict[b] for b in exp]) + else: + + arg_values = [evaluate_logic_exp(arg, bits_dict) for arg in exp.args] + + if exp.op in [BitWiseOp.AND, RegWiseOp.AND]: + return arg_values[0] & arg_values[1] + elif exp.op in [BitWiseOp.OR, RegWiseOp.OR]: + return arg_values[0] | arg_values[1] + elif exp.op in [BitWiseOp.XOR, RegWiseOp.XOR]: + return arg_values[0] ^ arg_values[1] + elif exp.op in [BitWiseOp.EQ, RegWiseOp.EQ]: + return int(arg_values[0] == arg_values[1]) + elif exp.op in [BitWiseOp.NEQ, RegWiseOp.NEQ]: + return int(arg_values[0] != arg_values[1]) + elif exp.op == BitWiseOp.NOT: + return 1 - arg_values[0] + elif exp.op == BitWiseOp.ZERO: + return 0 + elif exp.op == BitWiseOp.ONE: + return 1 + # elif exp.op == RegWiseOp.ADD: + # return arg_values[0] + arg_values[1] + # elif exp.op == RegWiseOp.SUB: + # return arg_values[0] - arg_values[1] + # elif exp.op == RegWiseOp.MUL: + # return arg_values[0] * arg_values[1] + # elif exp.op == RegWiseOp.POW: + # return int(arg_values[0] ** arg_values[1]) + # elif exp.op == RegWiseOp.LSH: + # return arg_values[0] << arg_values[1] + elif exp.op == RegWiseOp.RSH: + return arg_values[0] >> arg_values[1] + # elif exp.op == RegWiseOp.NEG: + # return -arg_values[0] + else: + # TODO: Currently not supporting RegWiseOp's DIV, EQ, NEQ, LT, GT, LEQ, + # GEQ and NOT, since these do not return int, so I am unsure what the + # semantic is meant to be. + # TODO: Similarly, it is not clear what to do with overflow of ADD, etc. + # so I have decided to not support them for now. + raise NotImplementedError( + f"Evaluation of {exp.op} not supported in ClassicalExpBox ", + "by pytket-cutensornet.", + ) + + +def from_little_endian(bitstring: list[bool]) -> int: + """Obtain the integer from the little-endian encoded bitstring (i.e. bitstring + [False, True] is interpreted as the integer 2).""" + return sum(1 << i for i, b in enumerate(bitstring) if b) diff --git a/pytket/extensions/cutensornet/structured_state/general.py b/pytket/extensions/cutensornet/structured_state/general.py index 1c8156af..c6d4276b 100644 --- a/pytket/extensions/cutensornet/structured_state/general.py +++ b/pytket/extensions/cutensornet/structured_state/general.py @@ -19,7 +19,14 @@ import numpy as np # type: ignore -from pytket.circuit import Command, Qubit +from pytket.circuit import ( + Command, + Op, + OpType, + Qubit, + Bit, + Conditional, +) from pytket.pauli import QubitPauliString try: @@ -28,6 +35,7 @@ warnings.warn("local settings failed to import cupy", ImportWarning) from pytket.extensions.cutensornet import CuTensorNetHandle +from .classical import apply_classical_command, from_little_endian # An alias for the CuPy type used for tensors try: @@ -47,7 +55,6 @@ def __init__( float_precision: Type[Any] = np.float64, value_of_zero: float = 1e-16, leaf_size: int = 8, - use_kahypar: bool = False, k: int = 4, optim_delta: float = 1e-5, loglevel: int = logging.WARNING, @@ -84,9 +91,6 @@ def __init__( ``np.float64`` precision (default) and ``1e-7`` for ``np.float32``. leaf_size: For ``TTN`` simulation only. Sets the maximum number of qubits in a leaf node when using ``TTN``. Default is 8. - use_kahypar: Use KaHyPar for graph partitioning (used in ``TTN``) if this - is True. Otherwise, use NetworkX (worse, but easy to setup). Defaults - to False. k: For ``MPSxMPO`` simulation only. Sets the maximum number of layers the MPO is allowed to have before being contracted. Increasing this might increase fidelity, but it will also increase resource requirements @@ -153,7 +157,6 @@ def __init__( raise ValueError("Maximum allowed leaf_size is 65.") self.leaf_size = leaf_size - self.use_kahypar = use_kahypar self.k = k self.optim_delta = 1e-5 self.loglevel = loglevel @@ -176,21 +179,19 @@ def copy(self) -> Config: class StructuredState(ABC): """Class representing a Tensor Network state.""" - @abstractmethod - def is_valid(self) -> bool: - """Verify that the tensor network state is valid. - - Returns: - False if a violation was detected or True otherwise. - """ - raise NotImplementedError(f"Method not implemented in {type(self).__name__}.") + _lib: CuTensorNetHandle + _cfg: Config + _logger: logging.Logger + _bits_dict: dict[Bit, bool] # Tracks the state of the classical variables - @abstractmethod def apply_gate(self, gate: Command) -> StructuredState: - """Applies the gate to the StructuredState. + """Apply the command to the `StructuredState`. + + Note: + Only one-qubit gates and two-qubit gates are supported. Args: - gate: The gate to be applied. + gate: The command to be applied. Returns: ``self``, to allow for method chaining. @@ -198,7 +199,65 @@ def apply_gate(self, gate: Command) -> StructuredState: Raises: RuntimeError: If the ``CuTensorNetHandle`` is out of scope. ValueError: If the command introduced is not a unitary gate. - ValueError: If gate acts on more than 2 qubits. + ValueError: If the command acts on more than 2 qubits. + """ + self._logger.debug(f"Applying {gate}.") + self._apply_command(gate.op, gate.qubits, gate.bits, gate.args) + return self + + def _apply_command( + self, op: Op, qubits: list[Qubit], bits: list[Bit], args: list[Any] + ) -> None: + """The implementation of `apply_gate`, acting on the unwrapped Command info.""" + if op.type == OpType.Measure: + q = qubits[0] + b = bits[0] + self._bits_dict[b] = self.measure({q}, destructive=False)[q] != 0 + + elif op.type == OpType.Reset: + assert len(qubits) + q = qubits[0] + # Measure and correct if outcome is |1> + outcome_1 = self.measure({q}, destructive=False)[q] != 0 + if outcome_1: + self._apply_command(Op.create(OpType.X), [q], [], [q]) + + elif op.is_gate(): # Either a unitary gate or a not supported "gate" + try: + unitary = op.get_unitary() + except: + raise ValueError(f"The command {op.type} introduced is not supported.") + + # Load the gate's unitary to the GPU memory + unitary = unitary.astype(dtype=self._cfg._complex_t, copy=False) + unitary = cp.asarray(unitary, dtype=self._cfg._complex_t) + + if len(qubits) not in [1, 2]: + raise ValueError( + "Gates must act on only 1 or 2 qubits! " + + f"This is not satisfied by {op.type}." + ) + + self.apply_unitary(unitary, qubits) + + elif isinstance(op, Conditional): + input_bits = args[: op.width] + tgt_value = op.value + # The input_bits encode a "value" int in little-endian + var_value = from_little_endian([self._bits_dict[b] for b in input_bits]) # type: ignore + # If the condition is apply the command in the body + if var_value == tgt_value: + self._apply_command(op.op, qubits, bits, args) + + else: # A purely classical operation + apply_classical_command(op, bits, args, self._bits_dict) + + @abstractmethod + def is_valid(self) -> bool: + """Verify that the tensor network state is valid. + + Returns: + False if a violation was detected or True otherwise. """ raise NotImplementedError(f"Method not implemented in {type(self).__name__}.") @@ -388,6 +447,10 @@ def get_amplitude(self, state: int) -> complex: """ raise NotImplementedError(f"Method not implemented in {type(self).__name__}.") + def get_bits(self) -> dict[Bit, bool]: + """Returns the dictionary of bits and their values.""" + return self._bits_dict.copy() + @abstractmethod def get_qubits(self) -> set[Qubit]: """Returns the set of qubits that ``self`` is defined on.""" diff --git a/pytket/extensions/cutensornet/structured_state/mps.py b/pytket/extensions/cutensornet/structured_state/mps.py index f62b5472..b9aee6e4 100644 --- a/pytket/extensions/cutensornet/structured_state/mps.py +++ b/pytket/extensions/cutensornet/structured_state/mps.py @@ -13,7 +13,7 @@ # limitations under the License. from __future__ import annotations # type: ignore import warnings -from typing import Union +from typing import Union, Optional from enum import Enum from random import Random # type: ignore @@ -29,7 +29,7 @@ except ImportError: warnings.warn("local settings failed to import cutensornet", ImportWarning) -from pytket.circuit import Command, Op, OpType, Qubit +from pytket.circuit import Op, OpType, Qubit, Bit from pytket.pauli import Pauli, QubitPauliString from pytket.extensions.cutensornet.general import CuTensorNetHandle, set_logger @@ -69,6 +69,7 @@ def __init__( libhandle: CuTensorNetHandle, qubits: list[Qubit], config: Config, + bits: Optional[list[Bit]] = None, ): """Initialise an MPS on the computational state ``|0>`` @@ -82,9 +83,6 @@ def __init__( tensor operations on the MPS. qubits: The list of qubits in the circuit to be simulated. config: The object describing the configuration for simulation. - - Raises: - ValueError: If less than two qubits are provided. """ self._lib = libhandle self._cfg = config @@ -93,11 +91,14 @@ def __init__( self._rng.seed(self._cfg.seed) self.fidelity = 1.0 + if bits is None: + self._bits_dict = dict() + else: + self._bits_dict = {b: False for b in bits} + n_tensors = len(qubits) if n_tensors == 0: # There's no initialisation to be done pass - elif n_tensors == 1: - raise ValueError("Please, provide at least two qubits.") else: self.qubit_position = {q: i for i, q in enumerate(qubits)} @@ -147,43 +148,6 @@ def is_valid(self) -> bool: return chi_ok and phys_ok and shape_ok and ds_ok - def apply_gate(self, gate: Command) -> MPS: - """Apply the gate to the MPS. - - Note: - Only one-qubit gates and two-qubit gates are supported. - - Args: - gate: The gate to be applied. - - Returns: - ``self``, to allow for method chaining. - - Raises: - RuntimeError: If the ``CuTensorNetHandle`` is out of scope. - ValueError: If the command introduced is not a unitary gate. - ValueError: If gate acts on more than 2 qubits. - """ - try: - unitary = gate.op.get_unitary() - except: - raise ValueError("The command introduced is not unitary.") - - # Load the gate's unitary to the GPU memory - unitary = unitary.astype(dtype=self._cfg._complex_t, copy=False) - unitary = cp.asarray(unitary, dtype=self._cfg._complex_t) - - self._logger.debug(f"Applying gate {gate}.") - if len(gate.qubits) not in [1, 2]: - raise ValueError( - "Gates must act on only 1 or 2 qubits! " - + f"This is not satisfied by {gate}." - ) - - self.apply_unitary(unitary, gate.qubits) - - return self - def apply_unitary( self, unitary: cp.ndarray, qubits: list[Qubit] ) -> StructuredState: @@ -269,6 +233,8 @@ def apply_qubit_relabelling(self, qubit_map: dict[Qubit, Qubit]) -> MPS: Raises: ValueError: If any of the keys in ``qubit_map`` are not qubits in the state. """ + self._flush() + new_qubit_position = dict() for q_orig, q_new in qubit_map.items(): # Check the qubit is in the state @@ -299,6 +265,8 @@ def add_qubit(self, new_qubit: Qubit, position: int, state: int = 0) -> MPS: ValueError: If ``position`` is negative or larger than ``len(self)``. ValueError: If ``state`` is not ``0`` or ``1``. """ + self._flush() + options = {"handle": self._lib.handle, "device_id": self._lib.device_id} if new_qubit in self.qubit_position.keys(): @@ -606,6 +574,7 @@ def measure(self, qubits: set[Qubit], destructive: bool = True) -> dict[Qubit, i Raises: ValueError: If an element in ``qubits`` is not a qubit in the state. """ + self._flush() result = dict() # Obtain the positions that need to be measured and build the reverse dict @@ -689,6 +658,8 @@ def postselect(self, qubit_outcomes: dict[Qubit, int]) -> float: ValueError: If all of the qubits in the MPS are being postselected. Instead, you may wish to use ``get_amplitude()``. """ + self._flush() + for q, v in qubit_outcomes.items(): if q not in self.qubit_position: raise ValueError(f"Qubit {q} is not a qubit in the MPS.") @@ -787,6 +758,8 @@ def expectation_value(self, pauli_string: QubitPauliString) -> float: Raises: ValueError: If a key in ``pauli_string`` is not a qubit in the MPS. """ + self._flush() + for q in pauli_string.map.keys(): if q not in self.qubit_position: raise ValueError(f"Qubit {q} is not a qubit in the MPS.") @@ -826,6 +799,7 @@ def expectation_value(self, pauli_string: QubitPauliString) -> float: def get_fidelity(self) -> float: """Returns the current fidelity of the state.""" + self._flush() return self.fidelity def get_statevector(self) -> np.ndarray: @@ -834,6 +808,8 @@ def get_statevector(self) -> np.ndarray: Raises: ValueError: If there are no qubits left in the MPS. """ + self._flush() + if len(self) == 0: raise ValueError("There are no qubits left in this MPS.") @@ -891,6 +867,7 @@ def get_amplitude(self, state: int) -> complex: Returns: The amplitude of the computational state in the MPS. """ + self._flush() # Auxiliar dictionary of physical bonds to qubit IDs qubit_id = {location: qubit for qubit, location in self.qubit_position.items()} diff --git a/pytket/extensions/cutensornet/structured_state/mps_mpo.py b/pytket/extensions/cutensornet/structured_state/mps_mpo.py index 3908b085..4dfbd9e5 100644 --- a/pytket/extensions/cutensornet/structured_state/mps_mpo.py +++ b/pytket/extensions/cutensornet/structured_state/mps_mpo.py @@ -28,7 +28,7 @@ except ImportError: warnings.warn("local settings failed to import cutensornet", ImportWarning) -from pytket.circuit import Qubit +from pytket.circuit import Qubit, Bit from pytket.extensions.cutensornet import CuTensorNetHandle from .general import Tensor, Config from .mps import ( @@ -49,6 +49,7 @@ def __init__( libhandle: CuTensorNetHandle, qubits: list[Qubit], config: Config, + bits: Optional[list[Bit]] = None, ): """Initialise an MPS on the computational state ``|0>``. @@ -63,7 +64,7 @@ def __init__( qubits: The list of qubits in the circuit to be simulated. config: The object describing the configuration for simulation. """ - super().__init__(libhandle, qubits, config) + super().__init__(libhandle, qubits, config, bits=bits) # Initialise the MPO data structure. This will keep a list of the gates # batched for application to the MPS; all of them will be applied at @@ -82,7 +83,7 @@ def __init__( # Initialise the MPS that we will use as first approximation of the # variational algorithm. - self._aux_mps = MPSxGate(libhandle, qubits, config) + self._aux_mps = MPSxGate(libhandle, qubits, config, bits=bits) self._mpo_bond_counter = 0 @@ -100,6 +101,47 @@ def update_libhandle(self, libhandle: CuTensorNetHandle) -> None: super().update_libhandle(libhandle) self._aux_mps.update_libhandle(libhandle) + def apply_qubit_relabelling(self, qubit_map: dict[Qubit, Qubit]) -> MPSxMPO: + """Relabels each qubit ``q`` as ``qubit_map[q]``. + + This does not apply any SWAP gate, nor it changes the internal structure of the + state. It simply changes the label of the physical bonds of the tensor network. + + Args: + qubit_map: Dictionary mapping each qubit to its new label. + + Returns: + ``self``, to allow for method chaining. + + Raises: + ValueError: If any of the keys in ``qubit_map`` are not qubits in the state. + """ + self._aux_mps.apply_qubit_relabelling(qubit_map) + super().apply_qubit_relabelling(qubit_map) + return self + + def add_qubit(self, new_qubit: Qubit, position: int, state: int = 0) -> MPS: + """Adds a qubit at the specified position. + + Args: + new_qubit: The identifier of the qubit to be added to the state. + position: The location the new qubit should be inserted at in the MPS. + Qubits on this and later indexed have their position shifted by 1. + state: Choose either ``0`` or ``1`` for the new qubit's state. + Defaults to ``0``. + + Returns: + ``self``, to allow for method chaining. + + Raises: + ValueError: If ``new_qubit`` already exists in the state. + ValueError: If ``position`` is negative or larger than ``len(self)``. + ValueError: If ``state`` is not ``0`` or ``1``. + """ + raise NotImplementedError( + "Creating new qubits is not currently supported in MPSxMPO." + ) + def _apply_1q_unitary(self, unitary: cp.ndarray, qubit: Qubit) -> MPSxMPO: """Applies the 1-qubit unitary to the MPS. @@ -117,6 +159,12 @@ def _apply_1q_unitary(self, unitary: cp.ndarray, qubit: Qubit) -> MPSxMPO: # Apply the gate to the MPS with eager approximation self._aux_mps._apply_1q_unitary(unitary, qubit) + if len(self) == 1: + # If there is only one qubit, there is no benefit in using an MPO, so + # simply copy from the standard MPS + self.tensors[0] = self._aux_mps.tensors[0].copy() + return self + # Glossary of bond IDs # i -> input to the MPO tensor # o -> output of the MPO tensor @@ -336,6 +384,10 @@ def _flush(self) -> None: The method applies variational optimisation of the MPS until it converges. Based on https://arxiv.org/abs/2207.05612. """ + if self._mpo_bond_counter == 0: + # MPO is empty, there is nothing to flush + return None + self._logger.info("Applying variational optimisation.") self._logger.info(f"Fidelity before optimisation={self._aux_mps.fidelity}") diff --git a/pytket/extensions/cutensornet/structured_state/simulation.py b/pytket/extensions/cutensornet/structured_state/simulation.py index 4f061877..e4756085 100644 --- a/pytket/extensions/cutensornet/structured_state/simulation.py +++ b/pytket/extensions/cutensornet/structured_state/simulation.py @@ -11,7 +11,7 @@ # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. -from typing import Optional +from typing import Optional, Any import warnings from enum import Enum @@ -26,7 +26,7 @@ except ImportError: warnings.warn("local settings failed to import kahypar", ImportWarning) -from pytket.circuit import Circuit, Command, Qubit +from pytket.circuit import Circuit, Command, OpType, Qubit from pytket.transform import Transform from pytket.architecture import Architecture from pytket.passes import DefaultMappingPass @@ -56,6 +56,7 @@ def simulate( circuit: Circuit, algorithm: SimulationAlgorithm, config: Config, + compilation_params: Optional[dict[str, Any]] = None, ) -> StructuredState: """Simulates the circuit and returns the ``StructuredState`` of the final state. @@ -73,6 +74,8 @@ def simulate( circuit: The pytket circuit to be simulated. algorithm: Choose between the values of the ``SimulationAlgorithm`` enum. config: The configuration object for simulation. + compilation_params: Experimental feature. Defaults to no compilation. + Parameters currently not documented. Returns: An instance of ``StructuredState`` for (an approximation of) the final state @@ -80,41 +83,59 @@ def simulate( """ logger = set_logger("Simulation", level=config.loglevel) - logger.info( - "Ordering the gates in the circuit to reduce canonicalisation overhead." - ) + if compilation_params is None: + compilation_params = dict() + + # Initialise the StructuredState if algorithm == SimulationAlgorithm.MPSxGate: state = MPSxGate( # type: ignore libhandle, circuit.qubits, - config, + bits=circuit.bits, + config=config, ) - sorted_gates = _get_sorted_gates(circuit, algorithm) elif algorithm == SimulationAlgorithm.MPSxMPO: state = MPSxMPO( # type: ignore libhandle, circuit.qubits, - config, + bits=circuit.bits, + config=config, ) - sorted_gates = _get_sorted_gates(circuit, algorithm) elif algorithm == SimulationAlgorithm.TTNxGate: + use_kahypar_option: bool = compilation_params.get("use_kahypar", False) + qubit_partition = _get_qubit_partition( - circuit, config.leaf_size, config.use_kahypar + circuit, config.leaf_size, use_kahypar_option ) state = TTNxGate( # type: ignore libhandle, qubit_partition, - config, + bits=circuit.bits, + config=config, ) - sorted_gates = _get_sorted_gates(circuit, algorithm, qubit_partition) + # If requested by the user, sort the gates to reduce canonicalisation overhead. + sort_gates_option: bool = compilation_params.get("sort_gates", False) + if sort_gates_option: + logger.info( + "Ordering the gates in the circuit to reduce canonicalisation overhead." + ) + + if algorithm == SimulationAlgorithm.TTNxGate: + commands = _get_sorted_gates(circuit, algorithm, qubit_partition) + else: + commands = _get_sorted_gates(circuit, algorithm) + else: + commands = circuit.get_commands() + + # Run the simulation logger.info("Running simulation...") # Apply the gates - for i, g in enumerate(sorted_gates): + for i, g in enumerate(commands): state.apply_gate(g) - logger.info(f"Progress... {(100*i) // len(sorted_gates)}%") + logger.info(f"Progress... {(100*i) // len(commands)}%") # Apply the batched operations that are left (if any) state._flush() @@ -152,6 +173,9 @@ def prepare_circuit_mps(circuit: Circuit) -> tuple[Circuit, dict[Qubit, Qubit]]: map of qubit names at the end of the circuit to their corresponding original names. """ + if circuit.n_qubits < 2: + # Nothing needs to be done + return (circuit, {q: q for q in circuit.qubits}) # Implement it in a line architecture cu = CompilationUnit(circuit) @@ -188,7 +212,7 @@ def _get_qubit_partition( A dictionary describing the partition in the format expected by TTN. Raises: - RuntimeError: If gate acts on more than 2 qubits. + RuntimeError: If a gate acts on more than 2 qubits. """ # Scan the circuit and generate the edges of the connectivity graph @@ -217,9 +241,9 @@ def _get_qubit_partition( # Apply balanced bisections until each qubit group is small enough partition = {0: circuit.qubits} - stop_bisec = False # Do at least one bisection (TTN reqs >1 leaf nodes) - while not stop_bisec: + # Stop if all groups have less than ``max_q_per_leaf`` qubits in them + while not all(len(group) <= max_q_per_leaf for group in partition.values()): old_partition = partition.copy() for key, group in old_partition.items(): # Apply the balanced bisection on this group @@ -235,9 +259,6 @@ def _get_qubit_partition( partition[2 * key] = groupA partition[2 * key + 1] = groupB - # Stop if all groups have less than ``max_q_per_leaf`` qubits in them - stop_bisec = all(len(group) <= max_q_per_leaf for group in partition.values()) - qubit_partition = {key: list(leaf_qubits) for key, leaf_qubits in partition.items()} return qubit_partition @@ -315,6 +336,13 @@ def _get_sorted_gates( 2-qubit gates that are close together are applied one after the other. This reduces the overhead of canonicalisation during simulation. + Notes: + If the circuit has any command (other than measurement) acting on bits, this + function gives up trying to sort the gates, and simply returns the standard + `circuit.get_commands()`. It would be possible to update this function so that + it can manage these commands as well, but it is not clear that there is a strong + use case for this. + Args: circuit: The original circuit. algorithm: The simulation algorithm that will be used on this circuit. @@ -326,6 +354,11 @@ def _get_sorted_gates( The same gates, ordered in a beneficial way for the given algorithm. """ all_gates = circuit.get_commands() + + # Abort if there is classical logic or classical control in the circuit (see note) + if any(len(g.bits) != 0 and g.op.type is not OpType.Measure for g in all_gates): + return all_gates + sorted_gates = [] # Entries from `all_gates` that are not yet in `sorted_gates` remaining = set(range(len(all_gates))) diff --git a/pytket/extensions/cutensornet/structured_state/ttn.py b/pytket/extensions/cutensornet/structured_state/ttn.py index 75dc531e..3a1a412e 100644 --- a/pytket/extensions/cutensornet/structured_state/ttn.py +++ b/pytket/extensions/cutensornet/structured_state/ttn.py @@ -30,7 +30,7 @@ except ImportError: warnings.warn("local settings failed to import cutensornet", ImportWarning) -from pytket.circuit import Command, Qubit +from pytket.circuit import Qubit, Bit from pytket.pauli import QubitPauliString from pytket.extensions.cutensornet.general import CuTensorNetHandle, set_logger @@ -96,6 +96,7 @@ def __init__( libhandle: CuTensorNetHandle, qubit_partition: dict[int, list[Qubit]], config: Config, + bits: Optional[list[Bit]] = None, ): """Initialise a TTN on the computational state ``|0>``. @@ -123,7 +124,6 @@ def __init__( ValueError: If the keys of ``qubit_partition`` do not range from ``0`` to ``2^l - 1`` for some ``l``. ValueError: If a ``Qubit`` is repeated in ``qubit_partition``. - ValueError: If there is only one entry in ``qubit_partition``. """ self._lib = libhandle self._cfg = config @@ -131,6 +131,11 @@ def __init__( self._rng = Random() self._rng.seed(self._cfg.seed) + if bits is None: + self._bits_dict = dict() + else: + self._bits_dict = {b: False for b in bits} + self.fidelity = 1.0 self.nodes: dict[RootPath, TreeNode] = dict() self.qubit_position: dict[Qubit, tuple[RootPath, int]] = dict() @@ -138,11 +143,6 @@ def __init__( n_groups = len(qubit_partition) if n_groups == 0: # There's no initialisation to be done pass - elif n_groups == 1: - raise ValueError( - "Only one entry to qubit_partition provided." - "Introduce a finer partition of qubits." - ) else: n_levels = math.floor(math.log2(n_groups)) if n_groups != 2**n_levels: @@ -242,37 +242,6 @@ def is_valid(self) -> bool: ) return chi_ok and phys_ok and rank_ok and shape_ok - def apply_gate(self, gate: Command) -> TTN: - """Apply the gate to the TTN. - - Note: - Only single-qubit gates and two-qubit gates are supported. - - Args: - gate: The gate to be applied. - - Returns: - ``self``, to allow for method chaining. - - Raises: - RuntimeError: If the ``CuTensorNetHandle`` is out of scope. - ValueError: If the command introduced is not a unitary gate. - ValueError: If gate acts on more than 2 qubits. - """ - try: - unitary = gate.op.get_unitary() - except: - raise ValueError("The command introduced is not unitary.") - - # Load the gate's unitary to the GPU memory - unitary = unitary.astype(dtype=self._cfg._complex_t, copy=False) - unitary = cp.asarray(unitary, dtype=self._cfg._complex_t) - - self._logger.debug(f"Applying gate {gate}.") - self.apply_unitary(unitary, gate.qubits) - - return self - def apply_unitary( self, unitary: cp.ndarray, qubits: list[Qubit] ) -> StructuredState: diff --git a/tests/conftest.py b/tests/conftest.py index a90aedf6..5297f495 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -51,6 +51,12 @@ def quantum_volume_circuit(n_qubits: int) -> Circuit: return c +@pytest.fixture +def q1_empty() -> Circuit: + circuit = Circuit(1) + return circuit + + @pytest.fixture def q5_empty() -> Circuit: circuit = Circuit(5) @@ -63,6 +69,13 @@ def q8_empty() -> Circuit: return circuit +@pytest.fixture +def q1_h0rz() -> Circuit: + circuit = Circuit(1) + circuit.H(0).Rz(0.3, 0) + return circuit + + @pytest.fixture def q2_x0() -> Circuit: circuit = Circuit(2) diff --git a/tests/test_structured_state.py b/tests/test_structured_state.py index 307db888..d3d78d61 100644 --- a/tests/test_structured_state.py +++ b/tests/test_structured_state.py @@ -248,8 +248,10 @@ def test_canonicalise_ttn(center: Union[RootPath, Qubit]) -> None: @pytest.mark.parametrize( "circuit", [ + pytest.lazy_fixture("q1_empty"), # type: ignore pytest.lazy_fixture("q5_empty"), # type: ignore pytest.lazy_fixture("q8_empty"), # type: ignore + pytest.lazy_fixture("q1_h0rz"), # type: ignore pytest.lazy_fixture("q2_x0"), # type: ignore pytest.lazy_fixture("q2_x1"), # type: ignore pytest.lazy_fixture("q2_v0"), # type: ignore @@ -308,7 +310,9 @@ def test_exact_circ_sim(circuit: Circuit, algorithm: SimulationAlgorithm) -> Non @pytest.mark.parametrize( "circuit", [ + pytest.lazy_fixture("q1_empty"), # type: ignore pytest.lazy_fixture("q5_empty"), # type: ignore + pytest.lazy_fixture("q1_h0rz"), # type: ignore pytest.lazy_fixture("q2_lcu1"), # type: ignore pytest.lazy_fixture("q2_lcu2"), # type: ignore pytest.lazy_fixture("q2_lcu3"), # type: ignore @@ -363,8 +367,10 @@ def test_prepare_circuit_mps(circuit: Circuit, algorithm: SimulationAlgorithm) - @pytest.mark.parametrize( "circuit", [ + pytest.lazy_fixture("q1_empty"), # type: ignore pytest.lazy_fixture("q5_empty"), # type: ignore pytest.lazy_fixture("q8_empty"), # type: ignore + pytest.lazy_fixture("q1_h0rz"), # type: ignore pytest.lazy_fixture("q2_x0"), # type: ignore pytest.lazy_fixture("q2_x1"), # type: ignore pytest.lazy_fixture("q2_v0"), # type: ignore @@ -409,8 +415,10 @@ def test_approx_circ_sim_gate_fid( @pytest.mark.parametrize( "circuit", [ + pytest.lazy_fixture("q1_empty"), # type: ignore pytest.lazy_fixture("q5_empty"), # type: ignore pytest.lazy_fixture("q8_empty"), # type: ignore + pytest.lazy_fixture("q1_h0rz"), # type: ignore pytest.lazy_fixture("q2_x0"), # type: ignore pytest.lazy_fixture("q2_x1"), # type: ignore pytest.lazy_fixture("q2_v0"), # type: ignore @@ -453,7 +461,9 @@ def test_approx_circ_sim_chi(circuit: Circuit, algorithm: SimulationAlgorithm) - @pytest.mark.parametrize( "circuit", [ + pytest.lazy_fixture("q1_empty"), # type: ignore pytest.lazy_fixture("q5_empty"), # type: ignore + pytest.lazy_fixture("q1_h0rz"), # type: ignore pytest.lazy_fixture("q2_x0cx01cx10"), # type: ignore pytest.lazy_fixture("q2_lcu2"), # type: ignore pytest.lazy_fixture("q3_cx01cz12x1rx0"), # type: ignore @@ -715,6 +725,7 @@ def test_expectation_value(circuit: Circuit, observable: QubitPauliString) -> No @pytest.mark.parametrize( "circuit", [ + pytest.lazy_fixture("q1_h0rz"), # type: ignore pytest.lazy_fixture("q2_v0cx01cx10"), # type: ignore pytest.lazy_fixture("q2_hadamard_test"), # type: ignore pytest.lazy_fixture("q2_lcu2"), # type: ignore diff --git a/tests/test_structured_state_conditionals.py b/tests/test_structured_state_conditionals.py new file mode 100644 index 00000000..c5b6a116 --- /dev/null +++ b/tests/test_structured_state_conditionals.py @@ -0,0 +1,505 @@ +import pytest + +import numpy as np + +from pytket.circuit import ( + Circuit, + CircBox, + OpType, + Qubit, + Bit, + if_not_bit, + reg_eq, +) +from pytket.circuit.logic_exp import BitWiseOp, create_bit_logic_exp + +from pytket.extensions.cutensornet.structured_state import ( + CuTensorNetHandle, + Config, + simulate, + SimulationAlgorithm, +) + + +# This first suite of tests comes from the pytket-qir extension +# (see https://github.com/CQCL/pytket-qir/blob/main/tests/conditional_test.py) +# Further down, there are tests to check that the simulation works correctly. + + +def test_circuit_with_classicalexpbox_i() -> None: + # test conditional handling + + circ = Circuit(3) + a = circ.add_c_register("a", 5) + b = circ.add_c_register("b", 5) + c = circ.add_c_register("c", 5) + d = circ.add_c_register("d", 5) + circ.H(0) + circ.add_classicalexpbox_register(a | b, c) # type: ignore + circ.add_classicalexpbox_register(c | b, d) # type: ignore + circ.add_classicalexpbox_register(c | b, d, condition=a[4]) # type: ignore + circ.H(0) + circ.Measure(Qubit(0), d[4]) + circ.H(1) + circ.Measure(Qubit(1), d[3]) + circ.H(2) + circ.Measure(Qubit(2), d[2]) + + with CuTensorNetHandle() as libhandle: + cfg = Config() + state = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, cfg) + assert state.is_valid() + assert np.isclose(state.vdot(state), 1.0, atol=cfg._atol) + assert state.get_fidelity() == 1.0 + + +def test_circuit_with_classicalexpbox_ii() -> None: + # test conditional handling with else case + + circ = Circuit(3) + a = circ.add_c_register("a", 5) + b = circ.add_c_register("b", 5) + c = circ.add_c_register("c", 5) + d = circ.add_c_register("d", 5) + circ.H(0) + circ.add_classicalexpbox_register(a | b, c) # type: ignore + circ.add_classicalexpbox_register(c | b, d) # type: ignore + circ.add_classicalexpbox_register( + c | b, d, condition=if_not_bit(a[4]) # type: ignore + ) + circ.H(0) + circ.Measure(Qubit(0), d[4]) + circ.H(1) + circ.Measure(Qubit(1), d[3]) + circ.H(2) + circ.Measure(Qubit(2), d[2]) + + with CuTensorNetHandle() as libhandle: + cfg = Config() + state = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, cfg) + assert state.is_valid() + assert np.isclose(state.vdot(state), 1.0, atol=cfg._atol) + assert state.get_fidelity() == 1.0 + + +@pytest.mark.skip(reason="Currently not supporting arithmetic operations in LogicExp") +def test_circuit_with_classicalexpbox_iii() -> None: + # test complicated conditions and recursive classical op + + circ = Circuit(2) + + a = circ.add_c_register("a", 15) + b = circ.add_c_register("b", 15) + c = circ.add_c_register("c", 15) + d = circ.add_c_register("d", 15) + e = circ.add_c_register("e", 15) + + circ.H(0) + bits = [Bit(i) for i in range(10)] + big_exp = bits[4] | bits[5] ^ bits[6] | bits[7] & bits[8] + circ.H(0, condition=big_exp) + + circ.add_classicalexpbox_register(a + b - d, c) # type: ignore + circ.add_classicalexpbox_register(a * b * d * c, e) # type: ignore + + with CuTensorNetHandle() as libhandle: + cfg = Config() + state = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, cfg) + assert state.is_valid() + assert np.isclose(state.vdot(state), 1.0, atol=cfg._atol) + assert state.get_fidelity() == 1.0 + + +def test_circuit_with_conditional_gate_i() -> None: + # test complicated conditions and recursive classical op + + circ = Circuit(2, 2).H(0).H(1).measure_all() + + circ.add_gate(OpType.H, [0], condition_bits=[0, 1], condition_value=3) + + with CuTensorNetHandle() as libhandle: + cfg = Config() + state = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, cfg) + assert state.is_valid() + assert np.isclose(state.vdot(state), 1.0, atol=cfg._atol) + assert state.get_fidelity() == 1.0 + + +def test_circuit_with_conditional_gate_ii() -> None: + # test complicated conditions and recursive classical op + + circ = Circuit(2, 3).H(0).H(1).measure_all() + + circ.add_gate(OpType.H, [0], condition_bits=[0, 1, 2], condition_value=3) + + with CuTensorNetHandle() as libhandle: + cfg = Config() + state = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, cfg) + assert state.is_valid() + assert np.isclose(state.vdot(state), 1.0, atol=cfg._atol) + assert state.get_fidelity() == 1.0 + + +def test_pcircuit_with_conditional_gate_iii() -> None: + # test conditional for manual added gates + + circ = Circuit(2, 3).H(0).H(1) + + circ.add_gate( + OpType.PhasedX, [0.1, 0.2], [0], condition_bits=[0, 1, 2], condition_value=3 + ) + + with CuTensorNetHandle() as libhandle: + cfg = Config() + state = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, cfg) + assert state.is_valid() + assert np.isclose(state.vdot(state), 1.0, atol=cfg._atol) + assert state.get_fidelity() == 1.0 + + +def test_circuit_with_conditional_gate_iv() -> None: + circ = Circuit(7, name="testcirc") + + syn = circ.add_c_register("syn", 4) + + circ.X(0, condition=reg_eq(syn, 1)) + circ.X(0, condition=reg_eq(syn, 2)) + circ.X(0, condition=reg_eq(syn, 2)) + circ.X(0, condition=reg_eq(syn, 3)) + circ.X(0, condition=reg_eq(syn, 4)) + circ.X(0, condition=reg_eq(syn, 4)) + + with CuTensorNetHandle() as libhandle: + cfg = Config() + state = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, cfg) + assert state.is_valid() + assert np.isclose(state.vdot(state), 1.0, atol=cfg._atol) + assert state.get_fidelity() == 1.0 + + +def test_pytket_qir_conditional_8() -> None: + c = Circuit(4) + c.H(0) + c.H(1) + c.H(2) + c.H(3) + cbox = CircBox(c) + d = Circuit(4) + a = d.add_c_register("a", 4) + d.add_circbox(cbox, [0, 2, 1, 3], condition=a[0]) + + with CuTensorNetHandle() as libhandle: + cfg = Config() + state = simulate(libhandle, c, SimulationAlgorithm.MPSxGate, cfg) + assert state.is_valid() + assert np.isclose(state.vdot(state), 1.0, atol=cfg._atol) + assert state.get_fidelity() == 1.0 + + +def test_pytket_qir_conditional_9() -> None: + c = Circuit(4) + c.X(0) + c.Y(1) + c.Z(2) + c.H(3) + cbox = CircBox(c) + d = Circuit(4) + a = d.add_c_register("a", 4) + d.add_circbox(cbox, [0, 2, 1, 3], condition=a[0]) + + with CuTensorNetHandle() as libhandle: + cfg = Config() + state = simulate(libhandle, c, SimulationAlgorithm.MPSxGate, cfg) + assert state.is_valid() + assert np.isclose(state.vdot(state), 1.0, atol=cfg._atol) + assert state.get_fidelity() == 1.0 + + +def test_pytket_qir_conditional_10() -> None: + box_circ = Circuit(4) + box_circ.X(0) + box_circ.Y(1) + box_circ.Z(2) + box_circ.H(3) + box_c = box_circ.add_c_register("c", 5) + + box_circ.H(0) + box_circ.add_classicalexpbox_register(box_c | box_c, box_c) # type: ignore + + cbox = CircBox(box_circ) + d = Circuit(4, 5) + a = d.add_c_register("a", 4) + d.add_circbox(cbox, [0, 2, 1, 3, 0, 1, 2, 3, 4], condition=a[0]) + + with CuTensorNetHandle() as libhandle: + cfg = Config() + state = simulate(libhandle, d, SimulationAlgorithm.MPSxGate, cfg) + assert state.is_valid() + assert np.isclose(state.vdot(state), 1.0, atol=cfg._atol) + assert state.get_fidelity() == 1.0 + + +def test_circuit_with_conditional_gate_v() -> None: + # test conditional with no register + + circ = Circuit(7, name="testcirc") + + exp = create_bit_logic_exp(BitWiseOp.ONE, []) + circ.H(0, condition=exp) + exp2 = create_bit_logic_exp(BitWiseOp.ZERO, []) + circ.H(0, condition=exp2) + + with CuTensorNetHandle() as libhandle: + cfg = Config() + state = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, cfg) + assert state.is_valid() + assert np.isclose(state.vdot(state), 1.0, atol=cfg._atol) + assert state.get_fidelity() == 1.0 + + +# The tests below check correctness of the simulator. + + +def test_correctness_reset_bits() -> None: + # This circuit does reset on two qubits. + n_shots = 10 + + circ = Circuit(2, 2).H(0).X(1).measure_all() + + circ.add_gate(OpType.X, [0], condition_bits=[0], condition_value=1) + circ.add_gate(OpType.X, [1], condition_bits=[1], condition_value=1) + + with CuTensorNetHandle() as libhandle: + cfg = Config() + + for _ in range(n_shots): + state = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, cfg) + assert state.is_valid() + assert np.isclose(state.vdot(state), 1.0, atol=cfg._atol) + assert state.get_fidelity() == 1.0 + # The outcome is the |00> state + assert np.isclose(abs(state.get_amplitude(0)), 1.0) + + +def test_correctness_reset_register() -> None: + # Test reset on register, including RangePredicate + n_shots = 10 + + circ = Circuit() + + q_reg = circ.add_q_register("a", 3) + circ.H(q_reg[0]) + circ.X(q_reg[1]) + circ.Rx(1.5, q_reg[2]) + + c_reg = circ.add_c_register("c", 3) + for q, c in zip(q_reg, c_reg): + circ.Measure(q, c) + + # Correct the least significant qubit (in an unnecessarily complicated way) + circ.add_gate(OpType.X, [q_reg[0]], condition_bits=c_reg, condition_value=1) + circ.add_gate(OpType.X, [q_reg[0]], condition_bits=c_reg, condition_value=3) + circ.add_gate(OpType.X, [q_reg[0]], condition_bits=c_reg, condition_value=5) + circ.add_gate(OpType.X, [q_reg[0]], condition_bits=c_reg, condition_value=7) + # Correct the middle qubit (straightforwad way) + circ.add_gate(OpType.X, [q_reg[1]], condition_bits=[c_reg[1]], condition_value=1) + # Correct the last bit using RangePredicate to create the flag + flag = circ.add_c_register("flag", 1) + circ.add_c_range_predicate( + minval=4, maxval=7, args_in=[b for b in c_reg], arg_out=flag[0] + ) + circ.add_gate(OpType.X, [q_reg[2]], condition_bits=flag, condition_value=1) + + with CuTensorNetHandle() as libhandle: + cfg = Config() + + for _ in range(n_shots): + state = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, cfg) + assert state.is_valid() + assert np.isclose(state.vdot(state), 1.0, atol=cfg._atol) + assert state.get_fidelity() == 1.0 + # The outcome is the |000> state + assert np.isclose(abs(state.get_amplitude(0)), 1.0) + + +def test_correctness_copy_bits() -> None: + # Dummy circuit where two bits are set and then copied + circ = Circuit(1) + orig = circ.add_c_register("orig", 2) + copied = circ.add_c_register("copied", 2) + # Set the `orig` bits to 0 and 1 + circ.add_c_setbits([False, True], [orig[0], orig[1]]) + # Copy the bits to the `copied` register + circ.add_c_copybits([orig[0], orig[1]], [copied[0], copied[1]]) + # Simulate + with CuTensorNetHandle() as libhandle: + cfg = Config() + state = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, cfg) + # Check that the copied register has the correct values + assert state.get_bits()[copied[0]] == False and state.get_bits()[copied[1]] == True + + +def test_correctness_teleportation_bit() -> None: + # A circuit to teleport a single qubit + + n_shots = 10 + + circ = Circuit(3, 2) + + # Generate an "interesting" state to be teleported + circ.Rx(0.42, 0) + + # Generate a Bell state + circ.H(1) + circ.CX(1, 2) + + # Apply Bell measurement + circ.CX(0, 1) + circ.H(0) + circ.Measure(0, 0) + circ.Measure(1, 1) + + # Apply conditional corrections + circ.add_gate(OpType.Z, [2], condition_bits=[0, 1], condition_value=1) + circ.add_gate(OpType.X, [2], condition_bits=[0, 1], condition_value=2) + circ.add_gate(OpType.Y, [2], condition_bits=[0, 1], condition_value=3) + + # Reset the other qubits + circ.add_gate(OpType.Reset, [0]) + circ.add_gate(OpType.Reset, [1]) + + with CuTensorNetHandle() as libhandle: + cfg = Config() + + for _ in range(n_shots): + state = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, cfg) + assert state.is_valid() + assert np.isclose(state.vdot(state), 1.0, atol=cfg._atol) + assert state.get_fidelity() == 1.0 + # The outcome is cos(0.42*pi/2) |000> - j*sin2(0.42*pi/2) |001> + assert np.isclose(abs(state.get_amplitude(0)) ** 2, 0.6243, atol=1e-4) + assert np.isclose(abs(state.get_amplitude(1)) ** 2, 0.3757, atol=1e-4) + + +def test_repeat_until_success_i() -> None: + # From Figure 8 of https://arxiv.org/pdf/1311.1074 + + attempts = 100 + + circ = Circuit() + qin = circ.add_q_register("qin", 1) + qaux = circ.add_q_register("aux", 1) + flag = circ.add_c_register("flag", 1) + circ.add_c_setbits([True], [flag[0]]) # Set flag bit to 1 + + for _ in range(attempts): + circ.add_gate(OpType.Reset, [qaux[0]], condition_bits=flag, condition_value=1) + circ.add_gate(OpType.H, [qaux[0]], condition_bits=flag, condition_value=1) + circ.add_gate(OpType.T, [qaux[0]], condition_bits=flag, condition_value=1) + circ.add_gate( + OpType.CX, [qaux[0], qin[0]], condition_bits=flag, condition_value=1 + ) + circ.add_gate(OpType.H, [qaux[0]], condition_bits=flag, condition_value=1) + circ.add_gate( + OpType.CX, [qaux[0], qin[0]], condition_bits=flag, condition_value=1 + ) + circ.add_gate(OpType.T, [qaux[0]], condition_bits=flag, condition_value=1) + circ.add_gate(OpType.H, [qaux[0]], condition_bits=flag, condition_value=1) + circ.Measure(qaux[0], flag[0], condition_bits=flag, condition_value=1) + + with CuTensorNetHandle() as libhandle: + cfg = Config() + + state = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, cfg) + assert state.is_valid() + assert np.isclose(state.vdot(state), 1.0, atol=cfg._atol) + assert state.get_fidelity() == 1.0 + + # The flag bit should have turned False + assert not state.get_bits()[flag[0]] + # The auxiliary qubits should be in state |0> + prob = state.postselect({qaux[0]: 0}) + assert np.isclose(prob, 1.0) + + target_state = [np.sqrt(1 / 3), np.sqrt(2 / 3) * 1j] + output_state = state.get_statevector() + # As indicated in the paper, the gate is implemented up to global phase + global_phase = target_state[0] / output_state[0] + assert np.isclose(abs(global_phase), 1.0) + output_state *= global_phase + assert np.allclose(target_state, output_state) + + +def test_repeat_until_success_ii() -> None: + # From Figure 1(c) of https://arxiv.org/pdf/1311.1074 + + attempts = 100 + + circ = Circuit() + qin = circ.add_q_register("qin", 1) + qaux = circ.add_q_register("aux", 2) + flag = circ.add_c_register("flag", 3) + circ.add_c_setbits([True, True], [flag[0], flag[1]]) # Set flag bits to 11 + circ.H(qin[0]) # Use to convert gate to sqrt(1/5)*I + i*sqrt(4/5)*X (i.e. Z -> X) + + for _ in range(attempts): + circ.add_classicalexpbox_bit( + flag[0] | flag[1], [flag[2]] + ) # Success if both are zero + + circ.add_gate( + OpType.Reset, [qaux[0]], condition_bits=[flag[2]], condition_value=1 + ) + circ.add_gate( + OpType.Reset, [qaux[1]], condition_bits=[flag[2]], condition_value=1 + ) + circ.add_gate(OpType.H, [qaux[0]], condition_bits=[flag[2]], condition_value=1) + circ.add_gate(OpType.H, [qaux[1]], condition_bits=[flag[2]], condition_value=1) + + circ.add_gate(OpType.T, [qin[0]], condition_bits=[flag[2]], condition_value=1) + circ.add_gate(OpType.Z, [qin[0]], condition_bits=[flag[2]], condition_value=1) + circ.add_gate( + OpType.Tdg, [qaux[0]], condition_bits=[flag[2]], condition_value=1 + ) + circ.add_gate( + OpType.CX, [qaux[1], qaux[0]], condition_bits=[flag[2]], condition_value=1 + ) + circ.add_gate(OpType.T, [qaux[0]], condition_bits=[flag[2]], condition_value=1) + circ.add_gate( + OpType.CX, [qin[0], qaux[1]], condition_bits=[flag[2]], condition_value=1 + ) + circ.add_gate(OpType.T, [qaux[1]], condition_bits=[flag[2]], condition_value=1) + + circ.add_gate(OpType.H, [qaux[0]], condition_bits=[flag[2]], condition_value=1) + circ.add_gate(OpType.H, [qaux[1]], condition_bits=[flag[2]], condition_value=1) + circ.Measure(qaux[0], flag[0], condition_bits=[flag[2]], condition_value=1) + circ.Measure(qaux[1], flag[1], condition_bits=[flag[2]], condition_value=1) + + # From chat with Silas and exploring the RUS as a block matrix, we have noticed + # that the circuit is missing an X correction when this condition is satisfied + circ.add_classicalexpbox_bit(flag[0] ^ flag[1], [flag[2]]) + circ.add_gate(OpType.Z, [qin[0]], condition_bits=[flag[2]], condition_value=1) + + circ.H(qin[0]) # Use to convert gate to sqrt(1/5)*I + i*sqrt(4/5)*X (i.e. Z -> X) + + with CuTensorNetHandle() as libhandle: + cfg = Config() + + state = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, cfg) + assert state.is_valid() + assert np.isclose(state.vdot(state), 1.0, atol=cfg._atol) + assert state.get_fidelity() == 1.0 + + # All of the flag bits should have turned False + assert all(not state.get_bits()[bit] for bit in flag) + # The auxiliary qubits should be in state |0> + prob = state.postselect({qaux[0]: 0, qaux[1]: 0}) + assert np.isclose(prob, 1.0) + + target_state = [np.sqrt(1 / 5), np.sqrt(4 / 5) * 1j] + output_state = state.get_statevector() + # As indicated in the paper, the gate is implemented up to global phase + global_phase = target_state[0] / output_state[0] + assert np.isclose(abs(global_phase), 1.0) + output_state *= global_phase + assert np.allclose(target_state, output_state)