diff --git a/.github/workflows/build-with-nix.yml b/.github/workflows/build-with-nix.yml index ac6ed09c..cc801114 100644 --- a/.github/workflows/build-with-nix.yml +++ b/.github/workflows/build-with-nix.yml @@ -39,3 +39,6 @@ jobs: - name: Test pytket-cutensornet # impure is necessary due to nixgl usage (system-dependent cuda) run: nix run .#tests --impure --accept-flake-config + - name: Test example notebooks + # impure is necessary due to nixgl usage (system-dependent cuda) + run: nix run .#example-tests --impure --accept-flake-config diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml index 807d038c..84df82a8 100644 --- a/.github/workflows/lint.yml +++ b/.github/workflows/lint.yml @@ -30,4 +30,4 @@ jobs: black --check . - name: Run pylint run: | - pylint --recursive=y --ignore=ttn_tutorial.py,mps_tutorial.py */ + pylint --recursive=y pytket/extensions/cutensornet/ diff --git a/docs/api.rst b/docs/api.rst index 1e91e2c8..48ed0e2c 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -1,11 +1,6 @@ API documentation ----------------- -.. autoclass:: pytket.extensions.cutensornet.CuTensorNetHandle - - .. automethod:: destroy - - .. toctree:: modules/general_state.rst modules/structured_state.rst diff --git a/docs/changelog.rst b/docs/changelog.rst index a1fc3c73..d2c71682 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -3,6 +3,15 @@ Changelog ~~~~~~~~~ +Unreleased +---------- + +* New API: ``GeneralBraOpKet`` for exact calculation of arbitrary ```` values. Can be used to calculate inner products, expectation values and arbitrary matrix elements. +* New feature: both ``GeneralState`` and ``GeneralBraOpKet`` admit circuits with parameterised gates. +* New feature: ``GeneralState`` has a new method ``get_amplitude`` to obtain the amplitude of computational basis states. +* New feature: ``GeneralState`` and ``CuTensorNetShotsBackend`` now support RNG seeds for sampling. +* Deprecated ``TensorNetwork`` object. It is still available for the sake of backwards compatibility, but it has been removed from doc pages. + 0.9.0 (October 2024) --------------------- diff --git a/docs/modules/general_state.rst b/docs/modules/general_state.rst index 0c552812..d9c219f5 100644 --- a/docs/modules/general_state.rst +++ b/docs/modules/general_state.rst @@ -5,29 +5,16 @@ General state (exact) simulation .. autoclass:: pytket.extensions.cutensornet.general_state.GeneralState() - .. automethod:: __init__ + .. automethod:: sample + .. automethod:: get_amplitude .. automethod:: get_statevector .. automethod:: expectation_value - .. automethod:: sample .. automethod:: destroy -cuQuantum `contract` API interface -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -.. autoclass:: pytket.extensions.cutensornet.general_state.TensorNetwork - -.. autoclass:: pytket.extensions.cutensornet.general_state.PauliOperatorTensorNetwork - -.. autoclass:: pytket.extensions.cutensornet.general_state.ExpectationValueTensorNetwork - -.. autofunction:: pytket.extensions.cutensornet.general_state.tk_to_tensor_network - -.. autofunction:: pytket.extensions.cutensornet.general_state.measure_qubits_state - -.. autofunction:: pytket.extensions.cutensornet.general_state.get_operator_expectation_value - -.. autofunction:: pytket.extensions.cutensornet.general_state.get_circuit_overlap +.. autoclass:: pytket.extensions.cutensornet.general_state.GeneralBraOpKet() + .. automethod:: contract + .. automethod:: destroy Pytket backend ~~~~~~~~~~~~~~ diff --git a/docs/modules/structured_state.rst b/docs/modules/structured_state.rst index 94cb320f..aeafb5df 100644 --- a/docs/modules/structured_state.rst +++ b/docs/modules/structured_state.rst @@ -3,6 +3,13 @@ Structured state evolution .. automodule:: pytket.extensions.cutensornet.structured_state +Library handle +~~~~~~~~~~~~~~ + +.. autoclass:: pytket.extensions.cutensornet.CuTensorNetHandle + + .. automethod:: destroy + Simulation ~~~~~~~~~~ diff --git a/examples/check-examples b/examples/check-examples index a9e947cc..47a470c1 100755 --- a/examples/check-examples +++ b/examples/check-examples @@ -9,7 +9,12 @@ do p2j -o -t ${name}-gen.ipynb python/${name}.py cmp ${name}.ipynb ${name}-gen.ipynb rm ${name}-gen.ipynb - # TODO, add this when GPU is added to CI - # Run script: - # python python/${name}.py + + # run tests are performed in nix, allowing + # us to manage the testing environment in a + # reproducible way. + # + # See /nix-support/pytket-cutensornet.nix, + # in the derivation called + # run-pytket-cutensornet-examples. done diff --git a/examples/ci-tested-notebooks.txt b/examples/ci-tested-notebooks.txt index 1cfc10e5..ee440a89 100644 --- a/examples/ci-tested-notebooks.txt +++ b/examples/ci-tested-notebooks.txt @@ -1,2 +1,3 @@ mps_tutorial -ttn_tutorial \ No newline at end of file +ttn_tutorial +general_state_tutorial \ No newline at end of file diff --git a/examples/general_state_tutorial.ipynb b/examples/general_state_tutorial.ipynb new file mode 100644 index 00000000..9a041f99 --- /dev/null +++ b/examples/general_state_tutorial.ipynb @@ -0,0 +1 @@ +{"cells": [{"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["import numpy as np\n", "from sympy import Symbol\n", "from scipy.stats import unitary_group # type: ignore\n", "from pytket.circuit import Circuit, OpType, Unitary2qBox, Qubit, Bit\n", "from pytket.passes import DecomposeBoxes\n", "from pytket.utils import QubitPauliOperator\n", "from pytket._tket.pauli import Pauli, QubitPauliString\n", "from pytket.circuit.display import render_circuit_jupyter"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["from pytket.extensions.cutensornet.general_state import (\n", " GeneralState,\n", " GeneralBraOpKet,\n", ")\n", "from pytket.extensions.cutensornet.backends import CuTensorNetShotsBackend"]}, {"cell_type": "markdown", "metadata": {}, "source": ["# Introduction
\n", "This notebook is a guide on how to use the features provided in the `general_state` submodule of pytket-cutensornet. This submodule is a thin wrapper of CuTensorNet's `NetworkState`, allowing users to convert pytket circuits into tensor networks and use CuTensorNet's contraction path optimisation algorithm.
\n", "All simulations realised with this submodule are *exact*. Once the pytket circuit has been converted to a tensor network, the computation has two steps:
\n", " 1. *Contraction path optimisation*. Attempts to find an order of contracting pairs of tensors in which the the total number of FLOPs is minimised. No operation on the tensor network occurs at this point. Runs on CPU.
\n", " 2. *Tensor network contraction*. Uses the ordering of contractions found in the previous step evaluate the tensor network. Runs on GPU.
\n", "
\n", "**Reference**: The original contraction path optimisation algorithm that NVIDIA implemented on CuTensorNet: https://arxiv.org/abs/2002.01935"]}, {"cell_type": "markdown", "metadata": {}, "source": ["# `GeneralState`
\n", "The class `GeneralState` is used to convert a circuit into a tensor network and query information from the final state. Let's walk through a simple example."]}, {"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)\n", "my_circ.measure_all()"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["render_circuit_jupyter(my_circ)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["The first step is to convert our pytket circuit into a tensor network. This is straightforward:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["tn_state = GeneralState(my_circ)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["The variable `tn_state` now holds a tensor network representation of `my_circ`.
\n", "**Note**: Circuits must not have mid-circuit measurements or classical logic. The measurements at the end of the circuit are stripped and only considered when calling `tn_state.sample(n_shots)`.
\n", "We can now query information from the state. For instance, let's calculate the probability of in the qubits 0 and 3 agreeing in their outcome."]}, {"cell_type": "markdown", "metadata": {}, "source": ["First, let's generate `|x>` computational basis states where `q[0]` and `q[3]` agree on their values. We can do this with some bitwise operators and list comprehension.
\n", "**Note**: Remember that pytket uses \"increasing lexicographic order\" (ILO) for qubits, so `q[0]` is the most significant bit."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["selected_states = [\n", " x\n", " for x in range(2**my_circ.n_qubits)\n", " if ( # Iterate over all possible states\n", " x & int(\"10000\", 2) == 0\n", " and x & int(\"00010\", 2) == 0 # both qubits are 0 or...\n", " or x & int(\"10000\", 2) != 0\n", " and x & int(\"00010\", 2) != 0 # both qubits are 1\n", " )\n", "]"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We can now query the amplitude of all of these states and calculate the probability by summing their squared absolute values."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["amplitudes = []\n", "for x in selected_states:\n", " amplitudes.append(tn_state.get_amplitude(x))\n", "probability = sum(abs(a) ** 2 for a in amplitudes)\n", "print(f\"Probability: {probability}\")"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Of course, calculating probabilities by considering the amplitudes of all relevant states is not efficient in general, since we may need to calculate a number of amplitudes that scales exponentially with the number of qubits. An alternative is to use expectation values. In particular, all of the states in `selected_states` are +1 eigenvectors of the `ZIIZI` observable and, hence, we can calculate the probability `p` by solving the equation ` = (+1)p + (-1)(1-p)` using the fact that `ZIIZI` only has +1 and -1 eigenvalues."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["string_ZIIZI = QubitPauliString(\n", " my_circ.qubits, [Pauli.Z, Pauli.I, Pauli.I, Pauli.Z, Pauli.I]\n", ")\n", "observable = QubitPauliOperator({string_ZIIZI: 1.0})\n", "expectation_val = tn_state.expectation_value(observable).real\n", "exp_probability = (expectation_val + 1) / 2\n", "assert np.isclose(probability, exp_probability, atol=0.0001)\n", "print(f\"Probability: {exp_probability}\")"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Alternatively, we can estimate the probability by sampling."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["n_shots = 100000\n", "outcomes = tn_state.sample(n_shots)\n", "hit_count = 0\n", "for bit_tuple, count in outcomes.get_counts().items():\n", " if bit_tuple[0] == bit_tuple[3]:\n", " hit_count += count\n", "samp_probability = hit_count / n_shots\n", "assert np.isclose(probability, samp_probability, atol=0.01)\n", "print(f\"Probability: {samp_probability}\")"]}, {"cell_type": "markdown", "metadata": {}, "source": ["When we finish doing computations with the `tn_state` we must destroy it to free GPU memory."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["tn_state.destroy()"]}, {"cell_type": "markdown", "metadata": {}, "source": ["To avoid forgetting this final step, we recommend users call `GeneralState` (and `GeneralBraOpKet`) as context managers:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["with GeneralState(my_circ) as my_state:\n", " expectation_val = my_state.expectation_value(observable)\n", "print(expectation_val)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Using this syntax, `my_state` is automatically destroyed when the code exists the `with ...` block."]}, {"cell_type": "markdown", "metadata": {}, "source": ["# Parameterised circuits
\n", "Circuits that only differ on the parameters of their gates have the same tensor network topology and, hence, we may use the same contraction path for all of them."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["a, b, c = Symbol(\"a\"), Symbol(\"b\"), Symbol(\"c\")\n", "param_circ1 = Circuit(5)\n", "param_circ1.Ry(a, 3).Ry(0.27, 4).CX(4, 3).Ry(b, 2).Ry(0.21, 3)\n", "param_circ1.Ry(0.12, 0).Ry(a, 1)\n", "param_circ1.add_gate(OpType.CnX, [0, 1, 4]).add_gate(OpType.CnX, [4, 1, 3])\n", "param_circ1.X(0).X(1).add_gate(OpType.CnY, [0, 1, 2]).add_gate(OpType.CnY, [0, 4, 3]).X(\n", " 0\n", ").X(1)\n", "param_circ1.Ry(-b, 0).Ry(-c, 1)\n", "render_circuit_jupyter(param_circ1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We can pass a parameterised circuit to `GeneralState`. The value of the parameters is provided when calling methods of `GeneralState`. The contraction path is automatically reused on different calls to the same method."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["n_circs = 5\n", "with GeneralState(param_circ1) as param_state:\n", " for i in range(n_circs):\n", " symbol_map = {s: np.random.random() for s in [a, b, c]}\n", " exp_val = param_state.expectation_value(observable, symbol_map=symbol_map)\n", " print(f\"Expectation value for circuit {i}: {exp_val.real}\")"]}, {"cell_type": "markdown", "metadata": {}, "source": ["# `GeneralBraOpKet`
\n", "The `GeneralBraOpKet` can be used to calculate any number that can be represented as the result of some `` where `|bra>` and `|ket>` are the final states of pytket circuits, and `op` is a `QubitPauliOperator`. The circuits for `|bra>` and `|ket>` need not be the same."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["x, y, z = Symbol(\"x\"), Symbol(\"y\"), Symbol(\"z\")\n", "param_circ2 = Circuit(5)\n", "param_circ2.H(0)\n", "param_circ2.S(1)\n", "param_circ2.Rz(x * z, 2)\n", "param_circ2.Ry(y + x, 3)\n", "param_circ2.TK1(x, y, z, 4)\n", "param_circ2.TK2(z - y, z - x, (x + y) * z, 1, 3)\n", "symbol_map = {a: 2.1, b: 1.3, c: 0.7, x: 3.0, y: 1.6, z: -8.3}"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We can calculate inner products by providing no `op`:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["with GeneralBraOpKet(bra=param_circ2, ket=param_circ1) as braket:\n", " inner_prod = braket.contract(symbol_map=symbol_map)\n", "with GeneralBraOpKet(bra=param_circ1, ket=param_circ2) as braket:\n", " inner_prod_conj = braket.contract(symbol_map=symbol_map)\n", "assert np.isclose(np.conj(inner_prod), inner_prod_conj)\n", "print(f\" = {inner_prod}\")\n", "print(f\" = {inner_prod_conj}\")"]}, {"cell_type": "markdown", "metadata": {}, "source": ["And we are not constrained to Hermitian operators:"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["string_XZIXX = QubitPauliString(\n", " param_circ2.qubits, [Pauli.X, Pauli.Z, Pauli.I, Pauli.X, Pauli.X]\n", ")\n", "string_IZZYX = QubitPauliString(\n", " param_circ2.qubits, [Pauli.I, Pauli.Z, Pauli.Z, Pauli.Y, Pauli.X]\n", ")\n", "string_ZIZXY = QubitPauliString(\n", " param_circ2.qubits, [Pauli.Z, Pauli.I, Pauli.Z, Pauli.X, Pauli.Y]\n", ")\n", "operator = QubitPauliOperator(\n", " {string_XZIXX: -1.38j, string_IZZYX: 2.36, string_ZIZXY: 0.42j + 0.3}\n", ")\n", "with GeneralBraOpKet(bra=param_circ2, ket=param_circ1) as braket:\n", " value = braket.contract(operator, symbol_map=symbol_map)\n", "print(value)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["# Backends
\n", "We provide a pytket `Backend` to obtain shots using `GeneralState`."]}, {"cell_type": "markdown", "metadata": {}, "source": ["Let's consider a more challenging circuit"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["def random_circuit(n_qubits: int, n_layers: int) -> Circuit:\n", " \"\"\"Random quantum volume circuit.\"\"\"\n", " c = Circuit(n_qubits, n_qubits)\n", " for _ in range(n_layers):\n", " qubits = np.random.permutation([i for i in range(n_qubits)])\n", " qubit_pairs = [[qubits[i], qubits[i + 1]] for i in range(0, n_qubits - 1, 2)]\n", " for pair in qubit_pairs:\n", " # Generate random 4x4 unitary matrix.\n", " SU4 = unitary_group.rvs(4) # random unitary in SU4\n", " SU4 = SU4 / (np.linalg.det(SU4) ** 0.25)\n", " SU4 = np.matrix(SU4)\n", " c.add_unitary2qbox(Unitary2qBox(SU4), *pair)\n", " DecomposeBoxes().apply(c)\n", " return c"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Let's measure only three of the qubits.
\n", "**Note**: The complexity of this simulation increases exponentially with the number of qubits measured. Other factors leading to intractability are circuit depth and qubit connectivity."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["n_shots = 1000\n", "quantum_vol_circ = random_circuit(n_qubits=40, n_layers=5)\n", "quantum_vol_circ.Measure(Qubit(0), Bit(0))\n", "quantum_vol_circ.Measure(Qubit(1), Bit(1))\n", "quantum_vol_circ.Measure(Qubit(2), Bit(2))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["The `CuTensorNetShotsBackend` is used in the same way as any other pytket `Backend`."]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["backend = CuTensorNetShotsBackend()\n", "compiled_circ = backend.get_compiled_circuit(quantum_vol_circ)\n", "results = backend.run_circuit(compiled_circ, n_shots=n_shots)\n", "print(results.get_counts())"]}], "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/general_state_tutorial.py b/examples/python/general_state_tutorial.py new file mode 100644 index 00000000..bef9b47d --- /dev/null +++ b/examples/python/general_state_tutorial.py @@ -0,0 +1,192 @@ +import numpy as np +from sympy import Symbol +from scipy.stats import unitary_group # type: ignore +from pytket.circuit import Circuit, OpType, Unitary2qBox, Qubit, Bit +from pytket.passes import DecomposeBoxes +from pytket.utils import QubitPauliOperator +from pytket._tket.pauli import Pauli, QubitPauliString +from pytket.circuit.display import render_circuit_jupyter + +from pytket.extensions.cutensornet.general_state import ( + GeneralState, + GeneralBraOpKet, +) +from pytket.extensions.cutensornet.backends import CuTensorNetShotsBackend + +# # Introduction +# This notebook is a guide on how to use the features provided in the `general_state` submodule of pytket-cutensornet. This submodule is a thin wrapper of CuTensorNet's `NetworkState`, allowing users to convert pytket circuits into tensor networks and use CuTensorNet's contraction path optimisation algorithm. +# All simulations realised with this submodule are *exact*. Once the pytket circuit has been converted to a tensor network, the computation has two steps: +# 1. *Contraction path optimisation*. Attempts to find an order of contracting pairs of tensors in which the the total number of FLOPs is minimised. No operation on the tensor network occurs at this point. Runs on CPU. +# 2. *Tensor network contraction*. Uses the ordering of contractions found in the previous step evaluate the tensor network. Runs on GPU. +# +# **Reference**: The original contraction path optimisation algorithm that NVIDIA implemented on CuTensorNet: https://arxiv.org/abs/2002.01935 + +# # `GeneralState` +# The class `GeneralState` is used to convert a circuit into a tensor network and query information from the final state. Let's walk through a simple example. + +my_circ = Circuit(5) +my_circ.CX(3, 4) +my_circ.H(2) +my_circ.CZ(0, 1) +my_circ.ZZPhase(0.1, 4, 3) +my_circ.TK2(0.3, 0.5, 0.7, 2, 1) +my_circ.Ry(0.2, 0) +my_circ.measure_all() + +render_circuit_jupyter(my_circ) + +# The first step is to convert our pytket circuit into a tensor network. This is straightforward: +tn_state = GeneralState(my_circ) + +# The variable `tn_state` now holds a tensor network representation of `my_circ`. +# **Note**: Circuits must not have mid-circuit measurements or classical logic. The measurements at the end of the circuit are stripped and only considered when calling `tn_state.sample(n_shots)`. +# We can now query information from the state. For instance, let's calculate the probability of in the qubits 0 and 3 agreeing in their outcome. + +# First, let's generate `|x>` computational basis states where `q[0]` and `q[3]` agree on their values. We can do this with some bitwise operators and list comprehension. +# **Note**: Remember that pytket uses "increasing lexicographic order" (ILO) for qubits, so `q[0]` is the most significant bit. +selected_states = [ + x + for x in range(2**my_circ.n_qubits) + if ( # Iterate over all possible states + x & int("10000", 2) == 0 + and x & int("00010", 2) == 0 # both qubits are 0 or... + or x & int("10000", 2) != 0 + and x & int("00010", 2) != 0 # both qubits are 1 + ) +] + +# We can now query the amplitude of all of these states and calculate the probability by summing their squared absolute values. +amplitudes = [] +for x in selected_states: + amplitudes.append(tn_state.get_amplitude(x)) +probability = sum(abs(a) ** 2 for a in amplitudes) +print(f"Probability: {probability}") + +# Of course, calculating probabilities by considering the amplitudes of all relevant states is not efficient in general, since we may need to calculate a number of amplitudes that scales exponentially with the number of qubits. An alternative is to use expectation values. In particular, all of the states in `selected_states` are +1 eigenvectors of the `ZIIZI` observable and, hence, we can calculate the probability `p` by solving the equation ` = (+1)p + (-1)(1-p)` using the fact that `ZIIZI` only has +1 and -1 eigenvalues. +string_ZIIZI = QubitPauliString( + my_circ.qubits, [Pauli.Z, Pauli.I, Pauli.I, Pauli.Z, Pauli.I] +) +observable = QubitPauliOperator({string_ZIIZI: 1.0}) +expectation_val = tn_state.expectation_value(observable).real +exp_probability = (expectation_val + 1) / 2 +assert np.isclose(probability, exp_probability, atol=0.0001) +print(f"Probability: {exp_probability}") + +# Alternatively, we can estimate the probability by sampling. +n_shots = 100000 +outcomes = tn_state.sample(n_shots) +hit_count = 0 +for bit_tuple, count in outcomes.get_counts().items(): + if bit_tuple[0] == bit_tuple[3]: + hit_count += count +samp_probability = hit_count / n_shots +assert np.isclose(probability, samp_probability, atol=0.01) +print(f"Probability: {samp_probability}") + +# When we finish doing computations with the `tn_state` we must destroy it to free GPU memory. +tn_state.destroy() + +# To avoid forgetting this final step, we recommend users call `GeneralState` (and `GeneralBraOpKet`) as context managers: +with GeneralState(my_circ) as my_state: + expectation_val = my_state.expectation_value(observable) +print(expectation_val) + +# Using this syntax, `my_state` is automatically destroyed when the code exists the `with ...` block. + +# # Parameterised circuits +# Circuits that only differ on the parameters of their gates have the same tensor network topology and, hence, we may use the same contraction path for all of them. +a, b, c = Symbol("a"), Symbol("b"), Symbol("c") +param_circ1 = Circuit(5) +param_circ1.Ry(a, 3).Ry(0.27, 4).CX(4, 3).Ry(b, 2).Ry(0.21, 3) +param_circ1.Ry(0.12, 0).Ry(a, 1) +param_circ1.add_gate(OpType.CnX, [0, 1, 4]).add_gate(OpType.CnX, [4, 1, 3]) +param_circ1.X(0).X(1).add_gate(OpType.CnY, [0, 1, 2]).add_gate(OpType.CnY, [0, 4, 3]).X( + 0 +).X(1) +param_circ1.Ry(-b, 0).Ry(-c, 1) +render_circuit_jupyter(param_circ1) + +# We can pass a parameterised circuit to `GeneralState`. The value of the parameters is provided when calling methods of `GeneralState`. The contraction path is automatically reused on different calls to the same method. +n_circs = 5 +with GeneralState(param_circ1) as param_state: + for i in range(n_circs): + symbol_map = {s: np.random.random() for s in [a, b, c]} + exp_val = param_state.expectation_value(observable, symbol_map=symbol_map) + print(f"Expectation value for circuit {i}: {exp_val.real}") + + +# # `GeneralBraOpKet` +# The `GeneralBraOpKet` can be used to calculate any number that can be represented as the result of some `` where `|bra>` and `|ket>` are the final states of pytket circuits, and `op` is a `QubitPauliOperator`. The circuits for `|bra>` and `|ket>` need not be the same. +x, y, z = Symbol("x"), Symbol("y"), Symbol("z") +param_circ2 = Circuit(5) +param_circ2.H(0) +param_circ2.S(1) +param_circ2.Rz(x * z, 2) +param_circ2.Ry(y + x, 3) +param_circ2.TK1(x, y, z, 4) +param_circ2.TK2(z - y, z - x, (x + y) * z, 1, 3) +symbol_map = {a: 2.1, b: 1.3, c: 0.7, x: 3.0, y: 1.6, z: -8.3} + +# We can calculate inner products by providing no `op`: +with GeneralBraOpKet(bra=param_circ2, ket=param_circ1) as braket: + inner_prod = braket.contract(symbol_map=symbol_map) +with GeneralBraOpKet(bra=param_circ1, ket=param_circ2) as braket: + inner_prod_conj = braket.contract(symbol_map=symbol_map) +assert np.isclose(np.conj(inner_prod), inner_prod_conj) +print(f" = {inner_prod}") +print(f" = {inner_prod_conj}") + +# And we are not constrained to Hermitian operators: +string_XZIXX = QubitPauliString( + param_circ2.qubits, [Pauli.X, Pauli.Z, Pauli.I, Pauli.X, Pauli.X] +) +string_IZZYX = QubitPauliString( + param_circ2.qubits, [Pauli.I, Pauli.Z, Pauli.Z, Pauli.Y, Pauli.X] +) +string_ZIZXY = QubitPauliString( + param_circ2.qubits, [Pauli.Z, Pauli.I, Pauli.Z, Pauli.X, Pauli.Y] +) +operator = QubitPauliOperator( + {string_XZIXX: -1.38j, string_IZZYX: 2.36, string_ZIZXY: 0.42j + 0.3} +) +with GeneralBraOpKet(bra=param_circ2, ket=param_circ1) as braket: + value = braket.contract(operator, symbol_map=symbol_map) +print(value) + +# # Backends +# We provide a pytket `Backend` to obtain shots using `GeneralState`. + + +# Let's consider a more challenging circuit +def random_circuit(n_qubits: int, n_layers: int) -> Circuit: + """Random quantum volume circuit.""" + c = Circuit(n_qubits, n_qubits) + + for _ in range(n_layers): + qubits = np.random.permutation([i for i in range(n_qubits)]) + qubit_pairs = [[qubits[i], qubits[i + 1]] for i in range(0, n_qubits - 1, 2)] + + for pair in qubit_pairs: + # Generate random 4x4 unitary matrix. + SU4 = unitary_group.rvs(4) # random unitary in SU4 + SU4 = SU4 / (np.linalg.det(SU4) ** 0.25) + SU4 = np.matrix(SU4) + c.add_unitary2qbox(Unitary2qBox(SU4), *pair) + + DecomposeBoxes().apply(c) + return c + + +# Let's measure only three of the qubits. +# **Note**: The complexity of this simulation increases exponentially with the number of qubits measured. Other factors leading to intractability are circuit depth and qubit connectivity. +n_shots = 1000 +quantum_vol_circ = random_circuit(n_qubits=40, n_layers=5) +quantum_vol_circ.Measure(Qubit(0), Bit(0)) +quantum_vol_circ.Measure(Qubit(1), Bit(1)) +quantum_vol_circ.Measure(Qubit(2), Bit(2)) + +# The `CuTensorNetShotsBackend` is used in the same way as any other pytket `Backend`. +backend = CuTensorNetShotsBackend() +compiled_circ = backend.get_compiled_circuit(quantum_vol_circ) +results = backend.run_circuit(compiled_circ, n_shots=n_shots) +print(results.get_counts()) diff --git a/flake.nix b/flake.nix index 04b67de9..72eb5a54 100644 --- a/flake.nix +++ b/flake.nix @@ -43,9 +43,9 @@ in { packages = { default = pkgs.pytket-cutensornet; - cupy = pkgs.cupy'; pytket-cutensornet = pkgs.pytket-cutensornet; tests = pkgs.run-pytket-cutensornet-tests; + example-tests = pkgs.run-pytket-cutensornet-examples; }; devShells = { default = pkgs.mkShell { buildInputs = [ pkgs.pytket-cutensornet ]; }; diff --git a/nix-support/pytket-cutensornet.nix b/nix-support/pytket-cutensornet.nix index 4895c30d..c6398cd4 100644 --- a/nix-support/pytket-cutensornet.nix +++ b/nix-support/pytket-cutensornet.nix @@ -65,4 +65,33 @@ EOF export LD_LIBRARY_PATH; ${test-env}/bin/pytest -s ${../tests}; ''; + run-pytket-cutensornet-examples = let + example-env = super.python3.withPackages(ps: with ps; [ + self.pytket-cutensornet + matplotlib + numpy + networkx + ipython + nbmake + pytest + ]); + nixgl-bin = self.lib.getExe self.nixgl.auto.nixGLNvidia; + in super.writeShellScriptBin "run-pytket-cutensornet-examples" '' + HOME=$(mktemp -d); + export HOME; + NIXGL_PATH="$(${nixgl-bin} printenv LD_LIBRARY_PATH)"; + WSL_PATH="/usr/lib/wsl/lib"; + LD_LIBRARY_PATH="$NIXGL_PATH:$WSL_PATH:$LD_LIBRARY_PATH"; + export LD_LIBRARY_PATH; + + example_dir=${../examples}; + set -e; + for name in `cat ''${example_dir}/ci-tested-notebooks.txt`; + do + ${example-env}/bin/pytest \ + --nbmake \ + -p no:cacheprovider \ + $example_dir/$name.ipynb; + done; + ''; } diff --git a/pytket/extensions/cutensornet/backends/cutensornet_backend.py b/pytket/extensions/cutensornet/backends/cutensornet_backend.py index b09d33c0..397f51c3 100644 --- a/pytket/extensions/cutensornet/backends/cutensornet_backend.py +++ b/pytket/extensions/cutensornet/backends/cutensornet_backend.py @@ -15,7 +15,6 @@ """Methods to allow tket circuits to be run on the cuTensorNet simulator.""" from abc import abstractmethod -import warnings from typing import List, Union, Optional, Sequence from uuid import uuid4 @@ -24,7 +23,6 @@ from pytket.backends.backend import KwargTypes, Backend, BackendResult from pytket.backends.backendinfo import BackendInfo from pytket.backends.resulthandle import _ResultIdTuple -from pytket.extensions.cutensornet.general import CuTensorNetHandle from pytket.extensions.cutensornet.general_state import ( GeneralState, ) @@ -46,11 +44,6 @@ CustomPass, ) -try: - from cuquantum.cutensornet import StateAttribute, SamplerAttribute # type: ignore -except ImportError: - warnings.warn("local settings failed to import cuquantum", ImportWarning) - from .._metadata import __extension_version__, __extension_name__ @@ -203,40 +196,39 @@ def process_circuits( corresponding get_ method. Note: - Any element from the ``StateAttribute`` enum (see NVIDIA's CuTensorNet + Any element from the ``TNConfig`` enum (see NVIDIA's CuTensorNet API) can be provided as arguments to this method. For instance: - ``process_circuits(..., CONFIG_NUM_HYPER_SAMPLES=100)``. + ``process_circuits(..., tn_config={"num_hyper_samples": 100})``. Args: circuits: List of circuits to be submitted. n_shots: Number of shots in case of shot-based calculation. This should be ``None``, since this backend does not support shots. valid_check: Whether to check for circuit correctness. + tnconfig: Optional. A dict of cuTensorNet ``TNConfig`` keys and + their values. scratch_fraction: Optional. Fraction of free memory on GPU to allocate as - scratch space. Defaults to `0.75`. + scratch space; value between 0 and 1. Defaults to ``0.8``. Returns: Results handle objects. """ - scratch_fraction = float(kwargs.get("scratch_fraction", 0.75)) # type: ignore - attributes = { - k: v for k, v in kwargs.items() if k in StateAttribute._member_names_ - } + scratch_fraction = float(kwargs.get("scratch_fraction", 0.8)) # type: ignore + tnconfig = kwargs.get("tnconfig", dict()) # type: ignore circuit_list = list(circuits) if valid_check: self._check_all_circuits(circuit_list) handle_list = [] - with CuTensorNetHandle() as libhandle: - for circuit in circuit_list: - tn = GeneralState(circuit, libhandle) - sv = tn.get_statevector(attributes, scratch_fraction) - res_qubits = [qb for qb in sorted(circuit.qubits)] - handle = ResultHandle(str(uuid4())) - self._cache[handle] = { - "result": BackendResult(q_bits=res_qubits, state=sv) - } - handle_list.append(handle) + for circuit in circuit_list: + with GeneralState( + circuit, attributes=tnconfig, scratch_fraction=scratch_fraction + ) as tn: + sv = tn.get_statevector() + res_qubits = [qb for qb in sorted(circuit.qubits)] + handle = ResultHandle(str(uuid4())) + self._cache[handle] = {"result": BackendResult(q_bits=res_qubits, state=sv)} + handle_list.append(handle) return handle_list @@ -280,9 +272,9 @@ def process_circuits( corresponding get_ method. Note: - Any element from the ``SamplerAttribute`` enum (see NVIDIA's CuTensorNet + Any element from the ``TNConfig`` enum (see NVIDIA's CuTensorNet API) can be provided as arguments to this method. For instance: - ``process_circuits(..., CONFIG_NUM_HYPER_SAMPLES=100)``. + ``process_circuits(..., tn_config={"num_hyper_samples": 100})``. Args: circuits: List of circuits to be submitted. @@ -290,23 +282,19 @@ def process_circuits( Optionally, this can be a list of shots specifying the number of shots for each circuit separately. valid_check: Whether to check for circuit correctness. + seed: An optional RNG seed. Different calls to ``process_circuits`` with the + same seed will generate the same list of shot outcomes. + tnconfig: Optional. A dict of cuTensorNet ``TNConfig`` keys and + their values. scratch_fraction: Optional. Fraction of free memory on GPU to allocate as - scratch space. Defaults to `0.75`. + scratch space; value between 0 and 1. Defaults to ``0.8``. Returns: Results handle objects. """ - scratch_fraction = float(kwargs.get("scratch_fraction", 0.75)) # type: ignore - attributes = { - k: v for k, v in kwargs.items() if k in SamplerAttribute._member_names_ - } - - if "seed" in kwargs and kwargs["seed"] is not None: - # Current CuTensorNet does not support seeds for Sampler. I created - # a feature request in their repository. - raise NotImplementedError( - "The backend does not currently support user-defined seeds." - ) + scratch_fraction = float(kwargs.get("scratch_fraction", 0.8)) # type: ignore + tnconfig = kwargs.get("tnconfig", dict()) # type: ignore + seed = kwargs.get("seed", None) if n_shots is None: raise ValueError( @@ -321,14 +309,13 @@ def process_circuits( if valid_check: self._check_all_circuits(circuit_list) handle_list = [] - with CuTensorNetHandle() as libhandle: - for circuit, circ_shots in zip(circuit_list, all_shots): - tn = GeneralState(circuit, libhandle) - handle = ResultHandle(str(uuid4())) - self._cache[handle] = { - "result": tn.sample(circ_shots, attributes, scratch_fraction) - } - handle_list.append(handle) + for circuit, circ_shots in zip(circuit_list, all_shots): + handle = ResultHandle(str(uuid4())) + with GeneralState( + circuit, attributes=tnconfig, scratch_fraction=scratch_fraction + ) as tn: + self._cache[handle] = {"result": tn.sample(circ_shots, seed=seed)} + handle_list.append(handle) return handle_list diff --git a/pytket/extensions/cutensornet/general_state/__init__.py b/pytket/extensions/cutensornet/general_state/__init__.py index 04c3357b..bb5dadea 100644 --- a/pytket/extensions/cutensornet/general_state/__init__.py +++ b/pytket/extensions/cutensornet/general_state/__init__.py @@ -25,4 +25,4 @@ get_circuit_overlap, ) -from .tensor_network_state import GeneralState +from .tensor_network_state import GeneralState, GeneralBraOpKet diff --git a/pytket/extensions/cutensornet/general_state/tensor_network_state.py b/pytket/extensions/cutensornet/general_state/tensor_network_state.py index bedfc0f2..e89917d5 100644 --- a/pytket/extensions/cutensornet/general_state/tensor_network_state.py +++ b/pytket/extensions/cutensornet/general_state/tensor_network_state.py @@ -14,7 +14,7 @@ from __future__ import annotations import logging -from typing import Union, Optional, Tuple, Dict +from typing import Union, Optional, Any import warnings try: @@ -22,407 +22,249 @@ except ImportError: warnings.warn("local settings failed to import cupy", ImportWarning) import numpy as np -from sympy import Expr # type: ignore +from sympy import Expr, Symbol # type: ignore from numpy.typing import NDArray -from pytket.circuit import Circuit, Qubit, Bit, OpType -from pytket.extensions.cutensornet.general import CuTensorNetHandle, set_logger +from pytket.circuit import Circuit, Qubit, Bit, OpType, Op +from pytket.extensions.cutensornet.general import set_logger from pytket.utils import OutcomeArray from pytket.utils.operators import QubitPauliOperator from pytket.backends.backendresult import BackendResult try: - import cuquantum as cq # type: ignore - from cuquantum import cutensornet as cutn # type: ignore + from cuquantum.cutensornet.experimental import NetworkState, NetworkOperator # type: ignore except ImportError: warnings.warn("local settings failed to import cuquantum", ImportWarning) class GeneralState: - """Wrapper of cuTensorNet object for exact simulations via path optimisation.""" + """Wrapper of cuTensorNet's NetworkState for exact simulation of states. + + Constructs a tensor network for the output state of a pytket circuit. + The qubits are assumed to be initialised in the ``|0>`` state. + The object stores the *uncontracted* tensor network. + + Note: + Preferably used as ``with GeneralState(...) as state:`` so that GPU memory is + automatically released after execution. + + The ``circuit`` must not contain any ``CircBox`` or non-unitary command. + + Args: + circuit: A pytket circuit to be converted into a tensor network. + attributes: Optional. A dict of cuTensorNet ``TNConfig`` keys and + their values. + scratch_fraction: Optional. Fraction of free memory on GPU to allocate as + scratch space; value between 0 and 1. Defaults to ``0.8``. + loglevel: Internal logger output level. Use 30 for warnings only, 20 for + verbose and 10 for debug mode. + """ def __init__( self, circuit: Circuit, - libhandle: CuTensorNetHandle, + attributes: Optional[dict] = None, + scratch_fraction: float = 0.8, loglevel: int = logging.WARNING, ) -> None: - """Constructs a tensor network for the output state of a pytket circuit. - - The qubits are assumed to be initialised in the ``|0>`` state. - The resulting object stores the *uncontracted* tensor network. - - Note: - A ``libhandle`` is created via a ``with CuTensorNetHandle() as libhandle:`` - statement. The device where the ``GeneralState`` is stored will match the - one specified by the library handle. - - Note: - The ``circuit`` must not contain any ``CircBox`` or non-unitary command. - - Args: - circuit: A pytket circuit to be converted to a tensor network. - libhandle: An instance of a ``CuTensorNetHandle``. - loglevel: Internal logger output level. - """ self._logger = set_logger("GeneralState", loglevel) - self._lib = libhandle - libhandle.print_device_properties(self._logger) # Remove end-of-circuit measurements and keep track of them separately # It also resolves implicit swaps - self._circuit, self._measurements = _remove_meas_and_implicit_swaps(circuit) - # Identify each qubit with the index of the bond that represents it in the + circuit, self._measurements = _remove_meas_and_implicit_swaps(circuit) + # Identify each qubit with the index of the bond that represents it in the # tensor network stored in this GeneralState. Qubits are sorted in increasing # lexicographical order, which is the TKET standard. - self._qubit_idx_map = {q: i for i, q in enumerate(sorted(self._circuit.qubits))} + self._qubit_idx_map = {q: i for i, q in enumerate(sorted(circuit.qubits))} - num_qubits = self._circuit.n_qubits + self._phase = circuit.phase + self.n_qubits = circuit.n_qubits dim = 2 # We are always dealing with qubits, not qudits - qubits_dims = (dim,) * num_qubits # qubit size - self._logger.debug(f"Converting a quantum circuit with {num_qubits} qubits.") - data_type = cq.cudaDataType.CUDA_C_64F # for now let that be hard-coded - - self._state = cutn.create_state( - self._lib.handle, cutn.StatePurity.PURE, num_qubits, qubits_dims, data_type + qubits_dims = (dim,) * len(circuit.qubits) # qubit size + data_type = "complex128" # for now let that be hard-coded + + self.tn_state = NetworkState( + qubits_dims, + dtype=data_type, + config=attributes, + options={ + "memory_limit": f"{int(scratch_fraction*100)}%", + "logger": self._logger, + }, ) - self._gate_tensors = [] - # Append all gates to the TN - for com in self._circuit.get_commands(): - try: - gate_unitary = com.op.get_unitary() - except: - raise ValueError( - "All commands in the circuit must be unitary gates. The circuit " - f"contains {com}; no unitary matrix could be retrived for it." - ) - self._gate_tensors.append(_formatted_tensor(gate_unitary, com.op.n_qubits)) + # Maintain a dict of tensor_id->Op for symbolic Ops to be update when the user + # calls any evaluation function with the symbols specified. + self._symbolic_ops: dict[int, Op] = dict() + # Append all gates to the NetworkState + commands = circuit.get_commands() + for com in commands: + is_fixed = len(com.op.free_symbols()) == 0 + if is_fixed: + try: + gate_unitary = com.op.get_unitary() + except Exception as e: + raise ValueError( + "All commands in the circuit must be unitary gates. The circuit" + f" contains {com}; no unitary matrix could be retrived from it." + ) from e + else: + # Dummy unitary to be updated later with user specified paramaters + gate_unitary = np.identity(2**com.op.n_qubits, dtype="complex128") + gate_qubit_indices = tuple(self._qubit_idx_map[qb] for qb in com.qubits) + tensor_id = self.tn_state.apply_tensor_operator( + gate_qubit_indices, + _formatted_tensor(gate_unitary, com.op.n_qubits), + immutable=is_fixed, + unitary=True, + ) - cutn.state_apply_tensor_operator( - handle=self._lib.handle, - tensor_network_state=self._state, - num_state_modes=com.op.n_qubits, - state_modes=gate_qubit_indices, - tensor_data=self._gate_tensors[-1].data.ptr, - tensor_mode_strides=0, - immutable=1, - adjoint=0, - unitary=1, + if not is_fixed: + self._symbolic_ops[tensor_id] = com.op + + # If the circuit has no gates, apply one identity gate so that CuTensorNet does + # not panic due to no tensor operator in the NetworkState + if len(commands) == 0: + tensor_id = self.tn_state.apply_tensor_operator( + (0,), + _formatted_tensor(np.identity(2), 1), + immutable=True, + unitary=True, ) def get_statevector( self, - attributes: Optional[dict] = None, - scratch_fraction: float = 0.75, + symbol_map: Optional[dict[Symbol, float]] = None, on_host: bool = True, ) -> Union[cp.ndarray, np.ndarray]: """Contracts the circuit and returns the final statevector. Args: - attributes: Optional. A dict of cuTensorNet `StateAttribute` keys and - their values. - scratch_fraction: Optional. Fraction of free memory on GPU to allocate as - scratch space. + symbol_map: A dictionary where each element of the pytket circuit's + ``.free_symbols()`` is assigned a real number. on_host: Optional. If ``True``, converts cupy ``ndarray`` to numpy - ``ndarray``, copying it to host device (CPU). - Raises: - MemoryError: If there is insufficient workspace on GPU. + ``ndarray``, copying it to host device (CPU). Defaults to ``True``. Returns: Either a ``cupy.ndarray`` on a GPU, or a ``numpy.ndarray`` on a host device (CPU). Arrays are returned in a 1D shape. + + Raises: + ValueError: If not every free symbol in the circuit is assigned a + value in ``symbol_map``. """ + _update_tensors(self.tn_state, self._symbolic_ops, symbol_map) + + self._logger.debug("(Statevector) contracting the TN") + state_vector = self.tn_state.compute_state_vector() + sv = state_vector.flatten() # Convert to 1D + if on_host: + sv = cp.asnumpy(sv) + # Apply the phase from the circuit + sv *= np.exp(1j * np.pi * self._phase) + return sv + + def get_amplitude( + self, + state: int, + symbol_map: Optional[dict[Symbol, float]] = None, + ) -> complex: + """Returns the amplitude of the chosen computational state. - #################################### - # Configure the TN for contraction # - #################################### - if attributes is None: - attributes = dict() - attribute_pairs = [ - (getattr(cutn.StateAttribute, k), v) for k, v in attributes.items() - ] - - for attr, val in attribute_pairs: - attr_dtype = cutn.state_get_attribute_dtype(attr) - attr_arr = np.asarray(val, dtype=attr_dtype) - cutn.state_configure( - self._lib.handle, - self._state, - attr, - attr_arr.ctypes.data, - attr_arr.dtype.itemsize, - ) + Note: + The result is equivalent to ``state.get_statevector[b]``, but this method + is faster when querying a single amplitude (or just a few). - try: - ###################################### - # Allocate workspace for contraction # - ###################################### - stream = cp.cuda.Stream() - free_mem = self._lib.dev.mem_info[0] - scratch_size = int(scratch_fraction * free_mem) - scratch_space = cp.cuda.alloc(scratch_size) - self._logger.debug( - f"Allocated {scratch_size} bytes of scratch memory on GPU" - ) - work_desc = cutn.create_workspace_descriptor(self._lib.handle) - - cutn.state_prepare( - self._lib.handle, - self._state, - scratch_size, - work_desc, - stream.ptr, - ) - workspace_size_d = cutn.workspace_get_memory_size( - self._lib.handle, - work_desc, - cutn.WorksizePref.RECOMMENDED, - cutn.Memspace.DEVICE, - cutn.WorkspaceKind.SCRATCH, - ) + Args: + state: The integer whose bitstring describes the computational state. + The qubits in the bitstring are in increasing lexicographic order. + symbol_map: A dictionary where each element of the pytket circuit's + ``.free_symbols()`` is assigned a real number. - if workspace_size_d <= scratch_size: - cutn.workspace_set_memory( - self._lib.handle, - work_desc, - cutn.Memspace.DEVICE, - cutn.WorkspaceKind.SCRATCH, - scratch_space.ptr, - workspace_size_d, - ) - self._logger.debug( - f"Set {workspace_size_d} bytes of workspace memory out of the" - f" allocated scratch space." - ) + Returns: + The amplitude of the computational state in ``self``. - else: - raise MemoryError( - f"Insufficient workspace size on the GPU device {self._lib.dev.id}" - ) + Raises: + ValueError: If not every free symbol in the circuit is assigned a + value in ``symbol_map``. + """ + _update_tensors(self.tn_state, self._symbolic_ops, symbol_map) - ################### - # Contract the TN # - ################### - state_vector = cp.empty( - (2,) * self._circuit.n_qubits, dtype="complex128", order="F" - ) - cutn.state_compute( - self._lib.handle, - self._state, - work_desc, - (state_vector.data.ptr,), - stream.ptr, - ) - stream.synchronize() - sv = state_vector.flatten() - if on_host: - sv = cp.asnumpy(sv) - # Apply the phase from the circuit - sv *= np.exp(1j * np.pi * self._circuit.phase) - return sv - - finally: - cutn.destroy_workspace_descriptor(work_desc) # type: ignore - del scratch_space + self._logger.debug("(Amplitude) contracting the TN") + bitstring = format(state, "b").zfill(self.n_qubits) + amplitude = self.tn_state.compute_amplitude(bitstring) + # Apply the phase from the circuit + amplitude *= np.exp(1j * np.pi * self._phase) + return complex(amplitude) def expectation_value( self, operator: QubitPauliOperator, - attributes: Optional[dict] = None, - scratch_fraction: float = 0.75, + symbol_map: Optional[dict[Symbol, float]] = None, ) -> complex: """Calculates the expectation value of the given operator. Args: - operator: The operator whose expectation value is to be measured. - attributes: Optional. A dict of cuTensorNet `ExpectationAttribute` keys - and their values. - scratch_fraction: Optional. Fraction of free memory on GPU to allocate as - scratch space. - - Raises: - ValueError: If the operator acts on qubits not present in the circuit. + operator: The operator whose expectation value is to be calculated. + symbol_map: A dictionary where each element of the pytket circuit's + ``.free_symbols()`` is assigned a real number. Returns: The expectation value. - """ - ############################################ - # Generate the cuTensorNet operator object # - ############################################ - pauli_tensors = { - "X": _formatted_tensor(np.asarray([[0, 1], [1, 0]]), 1), - "Y": _formatted_tensor(np.asarray([[0, -1j], [1j, 0]]), 1), - "Z": _formatted_tensor(np.asarray([[1, 0], [0, -1]]), 1), - "I": _formatted_tensor(np.asarray([[1, 0], [0, 1]]), 1), - } - num_qubits = self._circuit.n_qubits - qubits_dims = (2,) * num_qubits - data_type = cq.cudaDataType.CUDA_C_64F + Raises: + ValueError: If the operator acts on qubits not present in the circuit. + ValueError: If not every free symbol in the circuit is assigned a + value in ``symbol_map``. + """ + _update_tensors(self.tn_state, self._symbolic_ops, symbol_map) - tn_operator = cutn.create_network_operator( - self._lib.handle, num_qubits, qubits_dims, data_type - ) + self._logger.debug("(Expectation value) converting operator to NetworkOperator") - self._logger.debug("Adding operator terms:") - for pauli_string, coeff in operator._dict.items(): - if isinstance(coeff, Expr): - numeric_coeff = complex(coeff.evalf()) # type: ignore - else: - numeric_coeff = complex(coeff) # type: ignore - self._logger.debug(f" {numeric_coeff}, {pauli_string}") + paulis = ["I", "X", "Y", "Z"] + pauli_strs = dict() + for pstr, coeff in operator._dict.items(): # Raise an error if the operator acts on qubits that are not in the circuit - if any(q not in self._circuit.qubits for q in pauli_string.map.keys()): + if any(q not in self._qubit_idx_map.keys() for q in pstr.map.keys()): raise ValueError( - f"The operator is acting on qubits {pauli_string.map.keys()}, " + f"The operator is acting on qubits {pstr.map.keys()}, " "but some of these are not present in the circuit, whose set of " - f"qubits is: {self._circuit.qubits}." + f"qubits is: {self._qubit_idx_map.keys()}." ) - # Obtain the tensors corresponding to this operator - qubit_pauli_map = { - q: pauli_tensors[pauli.name] for q, pauli in pauli_string.map.items() - } - - num_pauli = len(qubit_pauli_map) - num_modes = (1,) * num_pauli - state_modes = tuple( - (self._qubit_idx_map[qb],) for qb in qubit_pauli_map.keys() - ) - gate_data = tuple(tensor.data.ptr for tensor in qubit_pauli_map.values()) - - cutn.network_operator_append_product( - handle=self._lib.handle, - tensor_network_operator=tn_operator, - coefficient=numeric_coeff, - num_tensors=num_pauli, - num_state_modes=num_modes, - state_modes=state_modes, - tensor_mode_strides=0, - tensor_data=gate_data, - ) - - ###################################################### - # Configure the cuTensorNet expectation value object # - ###################################################### - expectation = cutn.create_expectation( - self._lib.handle, self._state, tn_operator - ) + pauli_list = [pstr[q] for q in self._qubit_idx_map.keys()] + this_pauli_string = "".join(map(lambda x: paulis[x], pauli_list)) + pauli_strs[this_pauli_string] = complex(coeff) - if attributes is None: - attributes = dict() - attribute_pairs = [ - (getattr(cutn.ExpectationAttribute, k), v) for k, v in attributes.items() - ] - - for attr, val in attribute_pairs: - attr_dtype = cutn.expectation_get_attribute_dtype(attr) - attr_arr = np.asarray(val, dtype=attr_dtype) - cutn.expectation_configure( - self._lib.handle, - expectation, - attr, - attr_arr.ctypes.data, - attr_arr.dtype.itemsize, - ) + tn_operator = NetworkOperator.from_pauli_strings(pauli_strs, dtype="complex128") - try: - ###################################### - # Allocate workspace for contraction # - ###################################### - stream = cp.cuda.Stream() - free_mem = self._lib.dev.mem_info[0] - scratch_size = int(scratch_fraction * free_mem) - scratch_space = cp.cuda.alloc(scratch_size) - - self._logger.debug( - f"Allocated {scratch_size} bytes of scratch memory on GPU" - ) - work_desc = cutn.create_workspace_descriptor(self._lib.handle) - cutn.expectation_prepare( - self._lib.handle, - expectation, - scratch_size, - work_desc, - stream.ptr, - ) - workspace_size_d = cutn.workspace_get_memory_size( - self._lib.handle, - work_desc, - cutn.WorksizePref.RECOMMENDED, - cutn.Memspace.DEVICE, - cutn.WorkspaceKind.SCRATCH, - ) - - if workspace_size_d <= scratch_size: - cutn.workspace_set_memory( - self._lib.handle, - work_desc, - cutn.Memspace.DEVICE, - cutn.WorkspaceKind.SCRATCH, - scratch_space.ptr, - workspace_size_d, - ) - self._logger.debug( - f"Set {workspace_size_d} bytes of workspace memory out of the" - f" allocated scratch space." - ) - else: - raise MemoryError( - f"Insufficient workspace size on the GPU device {self._lib.dev.id}" - ) - - ################################# - # Compute the expectation value # - ################################# - expectation_value = np.empty(1, dtype="complex128") - state_norm = np.empty(1, dtype="complex128") - cutn.expectation_compute( - self._lib.handle, - expectation, - work_desc, - expectation_value.ctypes.data, - state_norm.ctypes.data, - stream.ptr, - ) - stream.synchronize() - - # Note: we can also return `state_norm.item()`, but this should be 1 since - # we are always running unitary circuits - assert np.isclose(state_norm.item(), 1.0) - - return expectation_value.item() # type: ignore - - finally: - ##################################################### - # Destroy the Operator and ExpectationValue objects # - ##################################################### - cutn.destroy_workspace_descriptor(work_desc) # type: ignore - cutn.destroy_expectation(expectation) - cutn.destroy_network_operator(tn_operator) - del scratch_space + self._logger.debug("(Expectation value) contracting the TN") + return complex(self.tn_state.compute_expectation(tn_operator)) def sample( self, n_shots: int, - attributes: Optional[dict] = None, - scratch_fraction: float = 0.75, + symbol_map: Optional[dict[Symbol, float]] = None, + seed: Optional[int] = None, ) -> BackendResult: """Obtains samples from the measurements at the end of the circuit. Args: n_shots: The number of samples to obtain. - attributes: Optional. A dict of cuTensorNet `SamplerAttribute` keys and - their values. - scratch_fraction: Optional. Fraction of free memory on GPU to allocate as - scratch space. - Raises: - ValueError: If the circuit contains no measurements. - MemoryError: If there is insufficient workspace on GPU. + symbol_map: A dictionary where each element of the pytket circuit's + ``.free_symbols()`` is assigned a real number. + seed: An optional RNG seed. Different calls to ``sample`` with the same + seed will generate the same list of shot outcomes. + Returns: A pytket ``BackendResult`` with the data from the shots. + + Raises: + ValueError: If the circuit contains no measurements. + ValueError: If not every free symbol in the circuit is assigned a + value in ``symbol_map``. """ + _update_tensors(self.tn_state, self._symbolic_ops, symbol_map) num_measurements = len(self._measurements) if num_measurements == 0: @@ -435,138 +277,328 @@ def sample( qbit_list, cbit_list = zip(*self._measurements.items()) measured_modes = tuple(self._qubit_idx_map[qb] for qb in qbit_list) - ############################################ - # Configure the cuTensorNet sampler object # - ############################################ + self._logger.debug("(Sampling) contracting the TN") + if seed is not None: + seed = abs(seed) # Must be a positive integer + samples = self.tn_state.compute_sampling( + nshots=n_shots, + modes=measured_modes, + seed=seed, + ) + + # Convert the data in `samples` to an `OutcomeArray` using `from_readouts` + # which expects a 2D array `samples[SampleId][QubitId]` of 0s and 1s. + self._logger.debug("(Sampling) converting samples to pytket Backend") + readouts = np.empty(shape=(n_shots, num_measurements), dtype=int) + sample_id = 0 + for bitstring, count in samples.items(): + outcome = [int(b) for b in bitstring] - sampler = cutn.create_sampler( - handle=self._lib.handle, - tensor_network_state=self._state, - num_modes_to_sample=num_measurements, - modes_to_sample=measured_modes, + for _ in range(count): + readouts[sample_id] = outcome + sample_id += 1 + + shots = OutcomeArray.from_readouts(readouts) + + # We need to specify which bits correspond to which columns in the shots + # table. Since cuTensorNet promises that the ordering of outcomes is + # determined by the ordering we provided as `measured_modes`, which in + # turn corresponds to the ordering of qubits in `qbit_list`, the fact that + # `cbit_list` has the appropriate order in relation to `self._measurements` + # implies this defines the ordering of classical bits we intend. + return BackendResult(c_bits=cbit_list, shots=shots) + + def destroy(self) -> None: + """Destroy the tensor network and free up GPU memory. + + The preferred approach is to use a context manager as in + ``with GeneralState(...) as state:``. Otherwise, the user must release + memory explicitly by calling ``destroy()``. + """ + self._logger.debug("Freeing memory of GeneralState") + self.tn_state.free() + + def __enter__(self) -> GeneralState: + return self + + def __exit__(self, exc_type: Any, exc_value: Any, exc_tb: Any) -> None: + self.destroy() + + +class GeneralBraOpKet: + """Constructs a tensor network for ````. + + The qubits in ``ket`` and ``bra`` are assumed to be initialised in the ``|0>`` + state. The object stores the *uncontracted* tensor network. + + Note: + Preferably used as ``with GeneralBraOpKet(...) as braket:`` so that GPU memory + is automatically released after execution. + + The circuits must not contain any ``CircBox`` or non-unitary command. + + The operator is provided when ``contract`` is called. + + Args: + bra: A pytket circuit describing the |bra> state. + ket: A pytket circuit describing the |ket> state. + attributes: Optional. A dict of cuTensorNet ``TNConfig`` keys and + their values. + scratch_fraction: Optional. Fraction of free memory on GPU to allocate as + scratch space; value between 0 and 1. Defaults to ``0.8``. + loglevel: Internal logger output level. Use 30 for warnings only, 20 for + verbose and 10 for debug mode. + + Raises: + ValueError: If the circuits for ``ket`` or ``bra`` contain measurements. + ValueError: If the set of qubits of ``ket`` and ``bra`` do not match. + """ + + def __init__( + self, + bra: Circuit, + ket: Circuit, + attributes: Optional[dict] = None, + scratch_fraction: float = 0.8, + loglevel: int = logging.WARNING, + ) -> None: + self._logger = set_logger("GeneralBraOpKet", loglevel) + + # Check that the circuits have the same qubits + if set(ket.qubits) != set(bra.qubits): + raise ValueError( + "The circuits given to GeneralBraOpKet must act on the same qubits." + ) + # Remove end-of-circuit measurements and keep track of them separately + # It also resolves implicit swaps + ket, meas = _remove_meas_and_implicit_swaps(ket) + if len(meas) != 0: + raise ValueError( + "The circuits given to a GeneralBraOpKet cannot have measurements." + ) + bra, meas = _remove_meas_and_implicit_swaps(bra) + if len(meas) != 0: + raise ValueError( + "The circuits given to a GeneralBraOpKet cannot have measurements." + ) + # Identify each qubit with the index of the bond that represents it in the + # tensor network stored in this GeneralState. Qubits are sorted in increasing + # lexicographical order, which is the TKET standard. + self._qubit_idx_map = {q: i for i, q in enumerate(sorted(ket.qubits))} + self.n_qubits = ket.n_qubits + + dim = 2 # We are always dealing with qubits, not qudits + qubits_dims = (dim,) * self.n_qubits # qubit size + data_type = "complex128" # for now let that be hard-coded + + self.tn = NetworkState( + qubits_dims, + dtype=data_type, + config=attributes, + options={ + "memory_limit": f"{int(scratch_fraction*100)}%", + "logger": self._logger, + }, ) - if attributes is None: - attributes = dict() - attribute_pairs = [ - (getattr(cutn.SamplerAttribute, k), v) for k, v in attributes.items() - ] - - for attr, val in attribute_pairs: - attr_dtype = cutn.sampler_get_attribute_dtype(attr) - attr_arr = np.asarray(val, dtype=attr_dtype) - cutn.sampler_configure( - self._lib.handle, - sampler, - attr, - attr_arr.ctypes.data, - attr_arr.dtype.itemsize, + # Maintain a dict of tensor_id->Op for symbolic Ops to be update when the user + # calls any evaluation function with the symbols specified. + self._symbolic_ops: dict[int, Op] = dict() + # Apply all commands from the ket circuit + self._logger.debug("Converting the ket circuit to a NetworkState") + commands = ket.get_commands() + for com in commands: + is_fixed = len(com.op.free_symbols()) == 0 + if is_fixed: + try: + gate_unitary = com.op.get_unitary() + except Exception as e: + raise ValueError( + "All commands in the circuit must be unitary gates. The circuit" + f" contains {com}; no unitary matrix could be retrived from it." + ) from e + else: + # Dummy unitary to be updated later with user specified paramaters + gate_unitary = np.identity(2**com.op.n_qubits, dtype="complex128") + + gate_qubit_indices = tuple(self._qubit_idx_map[qb] for qb in com.qubits) + tensor_id = self.tn.apply_tensor_operator( + gate_qubit_indices, + _formatted_tensor(gate_unitary, com.op.n_qubits), + immutable=is_fixed, + unitary=True, ) - try: - ###################################### - # Allocate workspace for contraction # - ###################################### - stream = cp.cuda.Stream() - free_mem = self._lib.dev.mem_info[0] - scratch_size = int(scratch_fraction * free_mem) - scratch_space = cp.cuda.alloc(scratch_size) - - self._logger.debug( - f"Allocated {scratch_size} bytes of scratch memory on GPU" + if not is_fixed: + self._symbolic_ops[tensor_id] = com.op + + # If the circuit has no gates, apply one identity gate so that CuTensorNet does + # not panic due to no tensor operator in the NetworkState + if len(commands) == 0: + tensor_id = self.tn.apply_tensor_operator( + (0,), + _formatted_tensor(np.identity(2), 1), + immutable=True, + unitary=True, ) - work_desc = cutn.create_workspace_descriptor(self._lib.handle) - cutn.sampler_prepare( - self._lib.handle, - sampler, - scratch_size, - work_desc, - stream.ptr, + self.ket_phase = ket.phase + + # Create a placeholder Pauli identity operator, to be replaced when calling + # contract. + self._logger.debug("Creating a placeholder operator and appending it") + self._pauli_op_ids = [] + for mode in range(self.n_qubits): + tensor_id = self.tn.apply_tensor_operator( + (mode,), + _formatted_tensor(np.identity(2), 1), + immutable=False, + unitary=True, ) - workspace_size_d = cutn.workspace_get_memory_size( - self._lib.handle, - work_desc, - cutn.WorksizePref.RECOMMENDED, - cutn.Memspace.DEVICE, - cutn.WorkspaceKind.SCRATCH, + self._pauli_op_ids.append(tensor_id) + + # Apply all commands from the adjoint of the bra circuit + self._logger.debug("Applying the dagger of the bra circuit to the NetworkState") + commands = bra.dagger().get_commands() + for com in commands: + is_fixed = len(com.op.free_symbols()) == 0 + if is_fixed: + try: + gate_unitary = com.op.get_unitary() + except Exception as e: + raise ValueError( + "All commands in the circuit must be unitary gates. The circuit" + f" contains {com}; no unitary matrix could be retrived from it." + ) from e + else: + # Dummy unitary to be updated later with user specified paramaters + gate_unitary = np.identity(2**com.op.n_qubits, dtype="complex128") + + gate_qubit_indices = tuple(self._qubit_idx_map[qb] for qb in com.qubits) + tensor_id = self.tn.apply_tensor_operator( + gate_qubit_indices, + _formatted_tensor(gate_unitary, com.op.n_qubits), + immutable=is_fixed, + unitary=True, ) - if workspace_size_d <= scratch_size: - cutn.workspace_set_memory( - self._lib.handle, - work_desc, - cutn.Memspace.DEVICE, - cutn.WorkspaceKind.SCRATCH, - scratch_space.ptr, - workspace_size_d, - ) - self._logger.debug( - f"Set {workspace_size_d} bytes of workspace memory out of the" - f" allocated scratch space." + if not is_fixed: + self._symbolic_ops[tensor_id] = com.op + + # If the circuit has no gates, apply one identity gate so that CuTensorNet does + # not panic due to no tensor operator in the NetworkState + if len(commands) == 0: + tensor_id = self.tn.apply_tensor_operator( + (0,), + _formatted_tensor(np.identity(2), 1), + immutable=True, + unitary=True, + ) + self.bra_phase = bra.phase + + def contract( + self, + operator: Optional[QubitPauliOperator] = None, + symbol_map: Optional[dict[Symbol, float]] = None, + ) -> complex: + """Contract the tensor network to obtain the value of ````. + + Args: + operator: A pytket ``QubitPauliOperator`` describing the operator. If not + given, then the identity operator is used, so it computes inner product. + symbol_map: A dictionary where each element of both pytket circuits' + ``.free_symbols()`` is assigned a real number. + + Returns: + The value of ````. + + Raises: + ValueError: If ``operator`` acts on qubits that are not in the circuits. + """ + _update_tensors(self.tn, self._symbolic_ops, symbol_map) + + paulis = ["I", "X", "Y", "Z"] + pauli_matrix = { + "I": np.identity(2), + "X": Op.create(OpType.X).get_unitary(), + "Y": Op.create(OpType.Y).get_unitary(), + "Z": Op.create(OpType.Z).get_unitary(), + } + pauli_strs: dict[str, complex] = dict() + + # Some care has to be taken when handling QubitPauliOperators, since identity + # Paulis may be omitted from the dictionary. + if operator is None: + pauli_strs = {"".join("I" for _ in range(self.n_qubits)): complex(1.0)} + else: + for tk_pstr, coeff in operator._dict.items(): + # Raise an error if the operator acts on qubits missing from the circuit + if any(q not in self._qubit_idx_map.keys() for q in tk_pstr.map.keys()): + raise ValueError( + f"The operator is acting on qubits {tk_pstr.map.keys()}, some " + "of these are missing from the set of qubits present in the " + f"circuits: {self._qubit_idx_map.keys()}." + ) + pauli_list = [tk_pstr[q] for q in self._qubit_idx_map.keys()] + this_pauli_string = "".join(map(lambda x: paulis[x], pauli_list)) + pauli_strs[this_pauli_string] = complex(coeff) + + # Calculate the value by iterating over all components of the QubitPauliOperator + value = 0.0 + zero_bitstring = "".join("0" for _ in range(self.n_qubits)) + for pstr, coeff in pauli_strs.items(): + # Update the NetworkState with this Pauli + self._logger.debug(f"Updating the tensors of the Pauli operator {pstr}") + for mode, pauli in enumerate(pstr): + self.tn.update_tensor_operator( + self._pauli_op_ids[mode], + _formatted_tensor(pauli_matrix[pauli], 1), + unitary=True, ) + + if isinstance(coeff, Expr): + numeric_coeff = complex(coeff.evalf()) # type: ignore else: - raise MemoryError( - f"Insufficient workspace size on the GPU device {self._lib.dev.id}" - ) + numeric_coeff = complex(coeff) - ########################### - # Sample from the circuit # - ########################### - samples = np.empty((num_measurements, n_shots), dtype="int64", order="F") - cutn.sampler_sample( - self._lib.handle, - sampler, - n_shots, - work_desc, - samples.ctypes.data, - stream.ptr, - ) - stream.synchronize() - - # Convert the data in `samples` to an `OutcomeArray` - # `samples` is a 2D numpy array `samples[SampleId][QubitId]`, which is - # the transpose of what `OutcomeArray.from_readouts` expects - shots = OutcomeArray.from_readouts(samples.T) - # We need to specify which bits correspond to which columns in the shots - # table. Since cuTensorNet promises that the ordering of outcomes is - # determined by the ordering we provided as `measured_modes`, which in - # turn corresponds to the ordering of qubits in `qbit_list`, the fact that - # `cbit_list` has the appropriate order in relation to `self._measurements` - # determines this defines the ordering of classical bits we intend. - return BackendResult(c_bits=cbit_list, shots=shots) - - finally: - ############################## - # Destroy the Sampler object # - ############################## - cutn.destroy_workspace_descriptor(work_desc) # type: ignore - cutn.destroy_sampler(sampler) - del scratch_space + # Compute the amplitude of the |0> state. Since NetworkState holds the + # circuit bra.dagger()*operator*ket|0>, the value of the amplitude at <0| + # will be . + self._logger.debug(f"Computing the contribution of Pauli operator {pstr}") + value += numeric_coeff * self.tn.compute_amplitude(zero_bitstring) + + # Apply the phases from the circuits + value *= np.exp(1j * np.pi * self.ket_phase) + value *= np.exp(-1j * np.pi * self.bra_phase) + + return complex(value) def destroy(self) -> None: """Destroy the tensor network and free up GPU memory. - Note: - Users are required to call `destroy()` when done using a - `GeneralState` object. GPU memory deallocation is not - guaranteed otherwise. + The preferred approach is to use a context manager as in + ``with GeneralBraOpKet(...) as braket:``. Otherwise, the user must release + memory explicitly by calling ``destroy()``. """ - cutn.destroy_state(self._state) + self._logger.debug("Freeing memory of GeneralBraOpKet") + self.tn.free() + + def __enter__(self) -> GeneralBraOpKet: + return self + + def __exit__(self, exc_type: Any, exc_value: Any, exc_tb: Any) -> None: + self.destroy() def _formatted_tensor(matrix: NDArray, n_qubits: int) -> cp.ndarray: """Convert a matrix to the tensor format accepted by NVIDIA's API.""" - # Transpose is needed because of the way cuTN stores tensors. - # See https://github.com/NVIDIA/cuQuantum/discussions/124 - # #discussioncomment-8683146 for details. - cupy_matrix = cp.asarray(matrix).T.astype(dtype="complex128", order="F") + cupy_matrix = cp.asarray(matrix, order="C").astype(dtype="complex128") # We also need to reshape since a matrix only has 2 bonds, but for an # n-qubit gate we want 2^n bonds for input and another 2^n for output - return cupy_matrix.reshape([2] * (2 * n_qubits), order="F") + return cupy_matrix.reshape([2] * (2 * n_qubits)) -def _remove_meas_and_implicit_swaps(circ: Circuit) -> Tuple[Circuit, Dict[Qubit, Bit]]: +def _remove_meas_and_implicit_swaps(circ: Circuit) -> tuple[Circuit, dict[Qubit, Bit]]: """Convert a pytket Circuit to an equivalent circuit with no measurements or implicit swaps. The measurements are returned as a map between qubits and bits. @@ -599,3 +631,44 @@ def _remove_meas_and_implicit_swaps(circ: Circuit) -> Tuple[Circuit, Dict[Qubit, pure_circ.add_phase(circ.phase) return pure_circ, measure_map # type: ignore + + +def _update_tensors( + tn: NetworkState, + symbolic_ops: dict[int, Op], + symbol_map: Optional[dict[Symbol, float]], +) -> None: + """Updates the tensors with the specified values for symbols. + + Args: + tn: The NetworkState that we intend to update. + symbolic_ops: A dictionary mapping ``tensor_id`` to the parameterised Op. + symbol_map: A dictionary mapping symbols to real values. + + Raises: + ValueError: If not every free symbol in the circuits is assigned a + value in ``symbol_map``. + """ + if symbol_map is None: + symbol_map = dict() + + for tensor_id, op in symbolic_ops.items(): + subs_params = op.params.copy() + if any(symb not in symbol_map.keys() for symb in op.free_symbols()): + raise ValueError( + f"Missing values for some of the free symbols {op.free_symbols()}. " + f"Symbols given: {symbol_map}." + ) + + for symb in op.free_symbols(): + subs_params = [ + p.subs(symb, symbol_map[symb]) if isinstance(p, Expr) else p + for p in subs_params + ] + gate_unitary = Op.create(op.type, subs_params).get_unitary() + + tn.update_tensor_operator( + tensor_id, + _formatted_tensor(gate_unitary, op.n_qubits), + unitary=True, + ) diff --git a/tests/conftest.py b/tests/conftest.py index 5297f495..650fc8be 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,5 +1,6 @@ import pytest import numpy as np +from sympy import Symbol from scipy.stats import unitary_group # type: ignore from pytket.circuit import Circuit, OpType, Unitary2qBox, ToffoliBox, Qubit from pytket.passes import DecomposeBoxes, CnXPairwiseDecomposition @@ -180,6 +181,20 @@ def q4_lcu1() -> Circuit: return circuit +@pytest.fixture +def q4_lcu1_parameterised() -> Circuit: + a, b, c = Symbol("a"), Symbol("b"), Symbol("c") + circuit = Circuit(4) + circuit.Ry(a, 3).Ry(0.27, 2).CX(2, 3).Ry(b, 2).Ry(0.21, 3) + circuit.Ry(0.12, 0).Ry(a, 1) + circuit.add_gate(OpType.CnX, [0, 1, 2]).add_gate(OpType.CnX, [0, 1, 3]) + circuit.X(0).X(1).add_gate(OpType.CnY, [0, 1, 2]).add_gate(OpType.CnY, [0, 1, 3]).X( + 0 + ).X(1) + circuit.Ry(-b, 0).Ry(-c, 1) + return circuit + + @pytest.fixture def q4_multicontrols() -> Circuit: circ = Circuit(4) @@ -229,6 +244,19 @@ def q5_h0s1rz2ry3tk4tk13() -> Circuit: return circuit +@pytest.fixture +def q5_h0s1rz2ry3tk4tk13_parameterised() -> Circuit: + a, b, c = Symbol("a"), Symbol("b"), Symbol("c") + circuit = Circuit(5) + circuit.H(0) + circuit.S(1) + circuit.Rz(a * c, 2) + circuit.Ry(b + a, 3) + circuit.TK1(a, b, c, 4) + circuit.TK2(a - b, c - a, (a + b) * c, 1, 3) + return circuit + + @pytest.fixture def q8_x0h2v5z6() -> Circuit: circuit = Circuit(8) diff --git a/tests/test_cutensornet_backend.py b/tests/test_cutensornet_backend.py index 25c9d31c..e58cf098 100644 --- a/tests/test_cutensornet_backend.py +++ b/tests/test_cutensornet_backend.py @@ -34,6 +34,22 @@ def test_sampler_bell() -> None: assert np.isclose(res.get_counts()[(1, 1)] / n_shots, 0.5, atol=0.01) +def test_sampler_bell_seed() -> None: + n_shots = 1000 + c = Circuit(2, 2) + c.H(0) + c.CX(0, 1) + c.measure_all() + b = CuTensorNetShotsBackend() + c = b.get_compiled_circuit(c) + res1 = b.run_circuit(c, n_shots=n_shots, seed=1234) + res2 = b.run_circuit(c, n_shots=n_shots, seed=1234) + res3 = b.run_circuit(c, n_shots=n_shots, seed=4321) + print(type(res1.get_shots())) + assert np.all(res1.get_shots() == res2.get_shots()) + assert np.any(res1.get_shots() != res3.get_shots()) + + def test_config_options() -> None: c = Circuit(2, 2) c.H(0) @@ -46,7 +62,7 @@ def test_config_options() -> None: h = b1.process_circuit( c, scratch_fraction=0.3, - CONFIG_NUM_HYPER_SAMPLES=100, + tn_config={"num_hyper_samples": 100}, ) assert np.allclose( b1.get_result(h).get_state(), np.asarray([1, 0, 0, 1]) * 1 / np.sqrt(2) @@ -60,7 +76,7 @@ def test_config_options() -> None: c, n_shots=n_shots, scratch_fraction=0.3, - CONFIG_NUM_HYPER_SAMPLES=100, + tn_config={"num_hyper_samples": 100}, ) assert res.get_shots().shape == (n_shots, 2) assert np.isclose(res.get_counts()[(0, 0)] / n_shots, 0.5, atol=0.01) diff --git a/tests/test_general_state.py b/tests/test_general_state.py index f3aadf5b..77ea5998 100644 --- a/tests/test_general_state.py +++ b/tests/test_general_state.py @@ -1,13 +1,13 @@ import random import numpy as np import pytest +from sympy import Symbol from pytket.circuit import Circuit, ToffoliBox, Qubit, Bit from pytket.passes import DecomposeBoxes, CnXPairwiseDecomposition from pytket.transform import Transform from pytket.pauli import QubitPauliString, Pauli from pytket.utils.operators import QubitPauliOperator -from pytket.extensions.cutensornet.general_state import GeneralState -from pytket.extensions.cutensornet.structured_state import CuTensorNetHandle +from pytket.extensions.cutensornet.general_state import GeneralState, GeneralBraOpKet @pytest.mark.parametrize( @@ -38,27 +38,32 @@ pytest.lazy_fixture("q8_x0h2v5z6"), # type: ignore ], ) -def test_get_statevec(circuit: Circuit) -> None: - with CuTensorNetHandle() as libhandle: - state = GeneralState(circuit, libhandle) - sv = state.get_statevector() +def test_basic_circs_state(circuit: Circuit) -> None: + sv_pytket = circuit.get_statevector() - sv_pytket = circuit.get_statevector() - assert np.allclose(sv, sv_pytket, atol=1e-10) + op = QubitPauliOperator( + { + QubitPauliString({q: Pauli.I for q in circuit.qubits}): 1.0, + } + ) - op = QubitPauliOperator( - { - QubitPauliString({q: Pauli.I for q in circuit.qubits}): 1.0, - } - ) + with GeneralState(circuit) as state: + sv = state.get_statevector() + assert np.allclose(sv, sv_pytket, atol=1e-10) # Calculate the inner product as the expectation value # of the identity operator: = - state = GeneralState(circuit, libhandle) ovl = state.expectation_value(op) assert ovl == pytest.approx(1.0) - state.destroy() + # Check that all amplitudes agree + for i in range(len(sv)): + assert np.isclose(sv[i], state.get_amplitude(i)) + + # Calculate the inner product again, using GeneralBraOpKet + with GeneralBraOpKet(circuit, circuit) as braket: + ovl = braket.contract() + assert ovl == pytest.approx(1.0) def test_sv_toffoli_box_with_implicit_swaps() -> None: @@ -83,10 +88,8 @@ def test_sv_toffoli_box_with_implicit_swaps() -> None: Transform.OptimiseCliffords().apply(ket_circ) # Convert and contract - with CuTensorNetHandle() as libhandle: - state = GeneralState(ket_circ, libhandle) + with GeneralState(ket_circ) as state: ket_net_vector = state.get_statevector() - state.destroy() # Compare to pytket statevector ket_pytket_vector = ket_circ.get_statevector() @@ -119,26 +122,23 @@ def to_bool_tuple(n_qubits: int, x: int) -> tuple: CnXPairwiseDecomposition().apply(ket_circ) Transform.OptimiseCliffords().apply(ket_circ) - with CuTensorNetHandle() as libhandle: - state = GeneralState(ket_circ, libhandle) + with GeneralState(ket_circ) as state: ket_net_vector = state.get_statevector() - ket_pytket_vector = ket_circ.get_statevector() - assert np.allclose(ket_net_vector, ket_pytket_vector) + ket_pytket_vector = ket_circ.get_statevector() + assert np.allclose(ket_net_vector, ket_pytket_vector) - # Calculate the inner product as the expectation value - # of the identity operator: = - op = QubitPauliOperator( - { - QubitPauliString({q: Pauli.I for q in ket_circ.qubits}): 1.0, - } - ) + # Calculate the inner product as the expectation value + # of the identity operator: = + op = QubitPauliOperator( + { + QubitPauliString({q: Pauli.I for q in ket_circ.qubits}): 1.0, + } + ) - state = GeneralState(ket_circ, libhandle) + with GeneralState(ket_circ) as state: ovl = state.expectation_value(op) - assert ovl == pytest.approx(1.0) - - state.destroy() + assert ovl == pytest.approx(1.0) @pytest.mark.parametrize( @@ -197,12 +197,16 @@ def test_expectation_value(circuit: Circuit, observable: QubitPauliOperator) -> exp_val_tket = observable.state_expectation(circuit.get_statevector()) # Calculate using GeneralState - with CuTensorNetHandle() as libhandle: - state = GeneralState(circuit, libhandle) + with GeneralState(circuit) as state: exp_val = state.expectation_value(observable) assert np.isclose(exp_val, exp_val_tket) - state.destroy() + + # Calculate using GeneralBraOpKet + with GeneralBraOpKet(circuit, circuit) as braket: + exp_val = braket.contract(observable) + + assert np.isclose(exp_val, exp_val_tket) @pytest.mark.parametrize( @@ -253,8 +257,7 @@ def test_sampler(circuit: Circuit, measure_all: bool) -> None: circuit.Measure(q, Bit(i)) # Sample using our library - with CuTensorNetHandle() as libhandle: - state = GeneralState(circuit, libhandle) + with GeneralState(circuit) as state: results = state.sample(n_shots) # Verify distribution matches theoretical probabilities @@ -276,4 +279,46 @@ def test_sampler(circuit: Circuit, measure_all: bool) -> None: assert np.isclose(count / n_shots, prob, atol=0.01) - state.destroy() + +@pytest.mark.parametrize( + "circuit", + [ + pytest.lazy_fixture("q4_lcu1_parameterised"), # type: ignore + pytest.lazy_fixture("q5_h0s1rz2ry3tk4tk13_parameterised"), # type: ignore + ], +) +@pytest.mark.parametrize( + "symbol_map", + [ + {Symbol("a"): 0.3, Symbol("b"): 0.42, Symbol("c"): -0.13}, + {Symbol("a"): 5.3, Symbol("b"): 1.42, Symbol("c"): -0.07, Symbol("d"): 0.53}, + ], +) +def test_parameterised(circuit: Circuit, symbol_map: dict[Symbol, float]) -> None: + state = GeneralState(circuit) + sv = state.get_statevector(symbol_map) + + circuit.symbol_substitution(symbol_map) + sv_pytket = circuit.get_statevector() + assert np.allclose(sv, sv_pytket, atol=1e-10) + + op = QubitPauliOperator( + { + QubitPauliString({q: Pauli.I for q in circuit.qubits}): 1.0, + } + ) + + # Calculate the inner product as the expectation value + # of the identity operator: = + with GeneralState(circuit) as state: + ovl = state.expectation_value(op) + assert ovl == pytest.approx(1.0) + + # Check that all amplitudes agree + for i in range(len(sv)): + assert np.isclose(sv[i], state.get_amplitude(i)) + + # Calculate the inner product again, using GeneralBraOpKet + with GeneralBraOpKet(circuit, circuit) as braket: + ovl = braket.contract() + assert ovl == pytest.approx(1.0)