{ "cells": [ { "cell_type": "markdown", "id": "970a3ec8", "metadata": {}, "source": [ "# Qrisp CUDA-Q Interface\n", "\n", "Ever wanted to combine Qrisp's expressive, Pythonic quantum programming with [CUDA-Q](https://nvidia.github.io/cuda-quantum/)'s high-performance simulation and execution? You're in the right place.\n", "\n", "\n", "This tutorial shows you how to run Qrisp programs on CUDA-Q through the [CUDA-Q interface](../../reference/Jasp/CUDA-Q%20Interface.rst). After a simple Bell state, you'll teleport a qubit using mid-circuit measurements, estimate the energy of a quantum magnet with iterative phase estimation, and solve a combinatorial optimization problem with the hybrid variational algorithm QAOA. To finish, you'll simulate the time evolution of a spin chain using block encodings, which are also a building block of many other quantum algorithms, such as linear system solvers. Every example compiles to a native CUDA-Q kernel.\n", "\n", "The interface ships as the separate package `qrisp-cudaq`. Install it with `pip install \"qrisp[cudaq]\"`.\n", "\n", "Let's dive in!" ] }, { "cell_type": "markdown", "id": "def462b6", "metadata": {}, "source": [ "## Hello World: Bell State\n", "\n", "Let's start with the Hello World of quantum computing: the classic Bell state example. This demonstrates how to define a Qrisp kernel, **run** it with CUDA-Q, and see the results. Notice how Pythonic everything ends up being." ] }, { "cell_type": "code", "execution_count": 1, "id": "6821716e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Single run result: 0\n", "10 shots: [3, 0, 0, 0, 0, 0, 0, 3, 3, 0]\n" ] } ], "source": [ "import warnings\n", "warnings.filterwarnings(\"ignore\")\n", "\n", "import cudaq\n", "from qrisp import *\n", "from qrisp_cudaq import cudaq_kernel\n", "\n", "@cudaq_kernel\n", "def bell():\n", " qv = QuantumVariable(2)\n", " h(qv[0])\n", " cx(qv[0], qv[1])\n", " return measure(qv)\n", "\n", "# Run the kernel once (statevector simulation)\n", "result = bell()\n", "print(\"Single run result:\", result)\n", "\n", "# Run the kernel multiple times (shots)\n", "results = cudaq.run(bell, shots_count=10)\n", "print(\"10 shots:\", results)" ] }, { "cell_type": "markdown", "id": "fa552018", "metadata": {}, "source": [ "The `run` [execution mode in CUDA-Q](https://nvidia.github.io/cuda-quantum/latest/using/examples/executing_kernels.html) naturally fits Qrisp's programming model and supports versatile features like real-time hybrid control flow. However, because most current hardware backends only support sampling, Qrisp's `cudaq_kernel` decorator allows you to pass the keyword argument `execution_mode=\"sample\"` (the default is `\"run\"`).\n", "\n", "Note that `sample` kernel functions must return `void`. The runtime aggregates the measurement results across all shots, returning a `SampleResult` histogram of the raw measurement bitstrings." ] }, { "cell_type": "code", "execution_count": 2, "id": "f505d209", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "10 shots (sample): { 00:5 11:5 }\n", "\n" ] } ], "source": [ "@cudaq_kernel(execution_mode=\"sample\")\n", "def bell():\n", " qv = QuantumVariable(2)\n", " h(qv[0])\n", " cx(qv[0], qv[1])\n", " measure(qv)\n", "\n", "results = cudaq.sample(bell, shots_count=10)\n", "print(\"10 shots (sample):\", results)" ] }, { "cell_type": "markdown", "id": "54c63512", "metadata": {}, "source": [ "The results contain bitstrings `00` and `11` for the Bell state $\\ket{\\psi}=\\frac{1}{\\sqrt{2}}(\\ket{00}+\\ket{11})$." ] }, { "cell_type": "markdown", "id": "6ca7af55", "metadata": {}, "source": [ "## Complex Examples: Hybrid Quantum-Classical Logic\n", "\n", "Let's shift into a higher gear! In a hybrid program, measurement results don't just come back at the end: they feed into classical logic that decides what the quantum computer does next, while the program is still running. Qrisp compiles this logic into the CUDA-Q kernel itself, so it runs in real time during every shot. This requires the default `run` execution mode; `sample` kernels can't use measurement results this way.\n", "\n", "The following three examples build on each other:\n", "\n", "- **Quantum teleportation** uses mid-circuit measurements to decide which corrections to apply (feed-forward).\n", "- **A dice game** runs a loop whose number of iterations depends on what the measurements return.\n", "- **Iterative phase estimation** combines both, a loop of rounds with measurement-dependent corrections, to estimate the energy of a quantum magnet. It also shows how ordinary Python libraries like NumPy and NetworkX help set up a kernel, together with [Qrisp's operators module](https://www.qrisp.eu/reference/Operators/index.html) and [trotterization](https://www.qrisp.eu/reference/Operators/generated/qrisp.operators.qubit.QubitOperator.trotterization.html#qrisp.operators.qubit.QubitOperator.trotterization) for Hamiltonian simulation.\n", "\n", "Keep one distinction in mind: plain Python code in a kernel, such as NumPy and NetworkX calls or `for` loops and `if` statements on Python values, is evaluated once, when the kernel is compiled, and decides which quantum operations end up in the kernel. A Python `for` loop, for example, is unrolled: its body is written into the kernel once per iteration, and all of these operations then run in every shot. Qrisp's control-flow constructs, such as [with control(...)](https://qrisp.eu/reference/Jasp/Control%20Flow/ClControl.html) on a measurement result, [q_while_loop](https://qrisp.eu/reference/Jasp/Control%20Flow/Prefix%20Control.html#qrisp.jasp.q_while_loop) and [q_fori_loop](https://qrisp.eu/reference/Jasp/Control%20Flow/Prefix%20Control.html#qrisp.jasp.q_fori_loop), are instead compiled into the kernel as real control flow and evaluated during every shot, which is what allows them to depend on measurement results.\n", "\n" ] }, { "cell_type": "markdown", "id": "bc58a28f", "metadata": {}, "source": [ "### Mid-Circuit Measurement with Feed-Forward\n", "\n", "A very common hybrid pattern is to measure some qubits in the middle of a computation and use the results to decide which operations come next. The textbook example is **quantum teleportation**: Alice sends the state of a qubit to Bob using a Bell pair they share. She measures her two qubits, and depending on her two results, Bob applies an `x` and/or a `z` correction to his qubit. Only with these corrections does Bob's qubit end up in the state that Alice sent.\n", "\n", "To check that it works, Bob undoes the state preparation at the end. If the teleportation succeeded, his qubit is back in $\\ket{0}$." ] }, { "cell_type": "code", "execution_count": 3, "id": "35c7fd22", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(m1, m2, Bob) for 8 shots: [(True, False, False), (True, False, False), (True, False, False), (False, True, False), (True, False, False), (False, False, False), (True, True, False), (True, False, False)]\n" ] } ], "source": [ "@cudaq_kernel\n", "def teleportation():\n", " message = QuantumBool() # Alice's qubit, holding the state to send\n", " alice = QuantumBool() # Alice's half of the shared Bell pair\n", " bob = QuantumBool() # Bob's half of the shared Bell pair\n", "\n", " # The state to teleport: cos(θ/2)|0> + sin(θ/2)|1>\n", " theta = 1.2\n", " ry(theta, message)\n", "\n", " # Alice and Bob share a Bell pair.\n", " h(alice)\n", " cx(alice, bob)\n", "\n", " # Alice measures her two qubits in the Bell basis.\n", " cx(message, alice)\n", " h(message)\n", " m1 = measure(message)\n", " m2 = measure(alice)\n", "\n", " # Feed-forward: Bob corrects his qubit based on Alice's results.\n", " with control(m2):\n", " x(bob)\n", " with control(m1):\n", " z(bob)\n", "\n", " # Undo the preparation: if the teleportation worked, Bob's qubit is back in |0>.\n", " ry(-theta, bob)\n", " return m1, m2, measure(bob)\n", "\n", "print(\"(m1, m2, Bob) for 8 shots:\", cudaq.run(teleportation, shots_count=8))\n" ] }, { "cell_type": "markdown", "id": "2029280a", "metadata": {}, "source": [ "Alice's results `m1` and `m2` are random: each of the four combinations shows up in about a quarter of the shots. Bob's result, however, is always `False`, i.e. his qubit is measured in $\\ket{0}$. Thanks to the corrections, his qubit holds exactly the state that Alice sent, whatever she measured. Without the two `with control(...)` blocks, Bob would measure `True` in half of the shots.\n", "\n", "The same pattern, measuring mid-circuit and correcting based on the results, is at the heart of adaptive protocols and quantum error correction.\n" ] }, { "cell_type": "markdown", "id": "ed29d2ed", "metadata": {}, "source": [ "### While-Loop with Accumulated Measurement\n", "\n", "This example demonstrates [q_while_loop](https://qrisp.eu/reference/Jasp/Control%20Flow/Prefix%20Control.html#qrisp.jasp.q_while_loop): a loop whose number of iterations is decided at runtime by measurement results.\n", "\n", "We play a simple dice game with an eight-sided die built from three qubits: keep rolling and add up the points until a roll comes up lower than 4. Putting the three qubits into a uniform superposition and measuring them gives a number between 0 and 7, each with equal probability; adding 1 turns it into a die face from 1 to 8. How many rolls a shot takes, and therefore how many points it collects, is not known in advance: it depends on what the measurements return.\n", "\n", "The loop carries four values: the last roll, the running total of accumulated points, the number of rolls and the qubits. Each iteration resets the qubits, rolls the die, adds the roll to the total and counts it. The loop continues as long as the last roll is 4 or higher, and the low roll that ends the loop still counts toward the total. The kernel returns the total and the number of rolls (more on multiple return values below).\n" ] }, { "cell_type": "code", "execution_count": 4, "id": "d7bf2bdf", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(total, rolls) for 10 shots: [(23, 4), (3, 1), (34, 7), (9, 2), (21, 4), (65, 12), (6, 2), (48, 8), (3, 1), (6, 2)]\n" ] } ], "source": [ "from qrisp.jasp import q_while_loop\n", "\n", "@cudaq_kernel\n", "def roll_until_low():\n", " qv = QuantumVariable(3)\n", "\n", " def body_fun(val):\n", " roll, total, rolls, qv = val\n", " reset(qv)\n", " h(qv)\n", " roll = measure(qv) + 1\n", " return roll, total + roll, rolls + 1, qv\n", "\n", " def cond_fun(val):\n", " roll, total, rolls, qv = val\n", " return roll >= 4\n", "\n", " # Start with a high roll so that the loop body runs at least once.\n", " roll, total, rolls, qv = q_while_loop(cond_fun, body_fun, (8, 0, 0, qv))\n", "\n", " return total, rolls\n", "\n", "print(\"(total, rolls) for 10 shots:\", cudaq.run(roll_until_low, shots_count=10))\n" ] }, { "cell_type": "markdown", "id": "1e3a93f5", "metadata": {}, "source": [ "Each roll ends the game with probability 3/8, so about 38% of shots stop after the first roll, while now and then a shot keeps going for ten rolls or more. On average a shot takes 8/3 ≈ 2.7 rolls and collects 12 points. The 5/3 rolls before the last one each land between 4 and 8, worth 6 points on average, which adds up to 10 points. The final roll, a 1, 2 or 3, adds another 2 points on average. Both the number of iterations and the accumulated total are decided at runtime by the measurement results." ] }, { "cell_type": "markdown", "id": "0d21ac7e", "metadata": {}, "source": [ "### Iterative Phase Estimation: Energy of a Heisenberg Ring\n", "\n", "Our example estimates the energy of a quantum antiferromagnet: the **Heisenberg model** on a ring of four spins, $H = \\sum_{(i,j)} (X_iX_j + Y_iY_j + Z_iZ_j)$, where the sum runs over neighboring spins. We build the ring with NetworkX, turn it into a Hamiltonian with Qrisp's operators module, and use the time evolution $U = e^{-iHt}$ from [trotterization](https://www.qrisp.eu/reference/Operators/generated/qrisp.operators.qubit.QubitOperator.trotterization.html#qrisp.operators.qubit.QubitOperator.trotterization) in **iterative phase estimation**. A single ancilla qubit is reused in two rounds, each measuring one bit of the energy. Before the second measurement, the ancilla gets a phase correction that is computed at runtime from the first result. The correction removes the phase contribution of the bit that was already measured, which would otherwise make the second result random. This is real-time feed-forward: classical measurement results steering the rest of the quantum computation. The rounds run in a [q_fori_loop](https://qrisp.eu/reference/Jasp/Control%20Flow/Prefix%20Control.html#qrisp.jasp.q_fori_loop), which compiles to a genuine loop in the CUDA-Q kernel: the kernel stays the same size no matter how many bits of the energy we read." ] }, { "cell_type": "code", "execution_count": 5, "id": "a2950b4b", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Energies for 10 shots: [-4.0, -8.0, -8.0, -8.0, -8.0, 4.0, -4.0, -8.0, -8.0, 4.0]\n", "Energy frequencies: [(-8.0, 317), (-4.0, 504), (0.0, 3), (4.0, 176)]\n", "Mean energy: -3.848\n" ] } ], "source": [ "from collections import Counter\n", "import networkx as nx\n", "import numpy as np\n", "from qrisp.jasp import q_fori_loop\n", "from qrisp.operators import X, Y, Z\n", "\n", "@cudaq_kernel\n", "def heisenberg_energy():\n", " # A ring of 4 spins; neighboring spins interact.\n", " G = nx.cycle_graph(4)\n", " # Heisenberg Hamiltonian on the ring, and its time evolution U = exp(-iHt)\n", " # via second-order Trotterization.\n", " H = sum(X(i) * X(j) + Y(i) * Y(j) + Z(i) * Z(j) for i, j in G.edges)\n", " U = H.trotterization(order=2)\n", "\n", " # Classical antiferromagnet (Néel state): neighboring spins point opposite ways.\n", " spins = QuantumVariable(G.number_of_nodes())\n", " for i in G.nodes:\n", " if i % 2:\n", " x(spins[i])\n", "\n", " # Iterative phase estimation: each round measures one bit of the phase\n", " # that U imprints, starting with the least significant bit.\n", " t = np.pi / 8\n", " n_bits = 2\n", " anc = QuantumBool()\n", "\n", " def iqpe_iteration(i, phi):\n", " # Rounds run from the highest power of U down to U itself.\n", " k = n_bits - 1 - i\n", " h(anc)\n", " with control(anc):\n", " # U^(2^k): repeat the time evolution 2^k times.\n", " U(spins, t=t, steps=2, iter=2**k)\n", " # Feed-forward: remove the contribution of the bits measured so far.\n", " p(-np.pi * phi, anc)\n", " h(anc)\n", " phi = (measure(anc) + phi) / 2\n", " reset(anc)\n", " return phi\n", "\n", " phi = q_fori_loop(0, n_bits, iqpe_iteration, 0.0)\n", "\n", " # The phase is phi = -E*t/(2*pi) mod 1; convert it back to the energy E.\n", " return 16 * (phi > 0.5) - 16 * phi\n", "\n", "print(\"Energies for 10 shots:\", cudaq.run(heisenberg_energy, shots_count=10))\n", "\n", "shots = cudaq.run(heisenberg_energy, shots_count=1000)\n", "print(\"Energy frequencies:\", sorted(Counter(float(e) for e in shots).items()))\n", "print(\"Mean energy:\", np.mean(shots))\n" ] }, { "cell_type": "markdown", "id": "bde4a74b", "metadata": {}, "source": [ "The Néel state, with neighboring spins pointing in opposite directions, is the classical picture of an antiferromagnet. Its energy is $-4$: each of the four $Z_iZ_j$ terms contributes $-1$, while the $X_iX_j$ and $Y_iY_j$ terms average to zero. But it is not an eigenstate of $H$, so phase estimation projects it onto eigenstates and returns their energies: $-8$ in about a third of the shots, $-4$ in half of them and $+4$ in a sixth. These fractions are the overlaps of the Néel state with the eigenstates ($1/3$, $1/2$ and $1/6$), and their average is again $-4$.\n", "\n", "The ground-state energy $-8$ lies far below the classical value of $-4$. The true ground state is an entangled superposition of spin configurations, and no direct measurement of the spins could reveal its energy; phase estimation is what turns it into a measurable number.\n", "\n", "Since the $X_iX_j$, $Y_iY_j$ and $Z_iZ_j$ terms do not commute, Trotterization only approximates $e^{-iHt}$. This small error is behind the rare result $0$; more Trotter steps make it rarer.\n", "\n", "Why $t = \\pi/8$? Phase estimation reads the phase $\\varphi$ in $U\\ket{\\psi} = e^{2\\pi i\\varphi}\\ket{\\psi}$. With $U = e^{-iHt}$, an eigenstate with energy $E$ has $\\varphi = -Et/(2\\pi) \\bmod 1$, which for $t = \\pi/8$ becomes $\\varphi = -E/16 \\bmod 1$. The energies of the ring ($-8$, $-4$, $0$ and $4$) are multiples of 4, so every phase is a multiple of $1/4$ and two rounds read it exactly. Two energies would only give the same phase if they differed by $2\\pi/t = 16$, which is more than the width of the spectrum, so every energy gets its own phase. The last line of the kernel inverts this relation: $E = -16\\varphi$, shifted up by 16 for phases above $1/2$, which belong to the positive energy.\n" ] }, { "cell_type": "markdown", "id": "ff331ca9", "metadata": {}, "source": [ "## Kernel Basics: Returns, Parameters, and Arrays\n", "\n", "Before diving into larger algorithms, it is worth covering three practical building blocks that come up in almost every non-trivial kernel:\n", "\n", "- **Multiple return values**: a kernel can return several measurement results at once, packed into a tuple.\n", "- **Parametrized kernels**: scalar values (`int`, `float`, `bool`) and fixed-length arrays can be passed into a kernel as classical parameters, enabling runtime-configurable circuits without recompilation.\n", "- **Static vs dynamic arrays**: arrays created inside the kernel from `np.array` are resolved at compile time and support NumPy transformations; arrays passed in as `FixedShapeNDArray` parameters or created with `jnp.array` are dynamic and are primarily used for indexed element access inside the kernel.\n", "\n", "The examples below illustrate each of these patterns in a minimal, executable form." ] }, { "cell_type": "markdown", "id": "442ff965", "metadata": {}, "source": [ "### Multiple Returns\n", "\n", "Multiple returns are supported; they are returned as a single tuple:" ] }, { "cell_type": "code", "execution_count": 6, "id": "609ff8e5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[(6.0, 3.0), (5.0, 2.0), (4.0, 1.0), (6.0, 3.0), (4.0, 1.0)]\n" ] } ], "source": [ "@cudaq_kernel\n", "def main():\n", " a = QuantumFloat(3)\n", " b = QuantumFloat(2)\n", " a[:] = 3\n", " h(b)\n", " a += b\n", " return measure(a), measure(b)\n", "\n", "print(cudaq.run(main, shots_count=5))" ] }, { "cell_type": "markdown", "id": "6da52ca6", "metadata": {}, "source": [ "The variable `a` is initialized in the state $\\ket{3}$, and `b` is in state $\\frac{1}{2}(\\ket{0}+\\ket{1}+\\ket{2}+\\ket{3})$. Then `b` is added to `a`, resulting in the entangled state $\\ket{a}\\ket{b}=\\frac{1}{2}(\\ket{3}\\ket{0}+\\ket{4}\\ket{1}+\\ket{5}\\ket{2}+\\ket{6}\\ket{3})$. " ] }, { "cell_type": "markdown", "id": "12815b6a", "metadata": {}, "source": [ "### Parametrized Kernels\n", "\n", "Kernel functions decorated with `@cudaq_kernel` can accept classical parameters. The supported scalar parameter types are `int`, `float`, and `bool`. For array-valued parameters, use `FixedShapeNDArray[dtype, length]`, where `dtype` is one of those scalar types.\n", "\n", "The example below shows a simple scalar-parametrized kernel, a single-qubit rotation by a runtime angle:" ] }, { "cell_type": "code", "execution_count": 7, "id": "6bfe82ee", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Parametrized kernel (single run): 1.0\n" ] } ], "source": [ "@cudaq_kernel\n", "def parametrized_kernel(angle: float):\n", " qv = QuantumFloat(1)\n", " ry(angle, qv[0])\n", " return measure(qv)\n", "\n", "angle = 1.57\n", "result_static = parametrized_kernel(angle)\n", "print(\"Parametrized kernel (single run):\", result_static)" ] }, { "cell_type": "markdown", "id": "58b4ac47", "metadata": {}, "source": [ "#### Array Parameters: Static vs Dynamic Indexing\n", "\n", "Array parameters are annotated with `FixedShapeNDArray[dtype, length]`. There are two important access patterns:\n", "\n", "- **Static indexing**: the index is known at compile time (for example `angles[0]`).\n", "- **Dynamic indexing**: the index is computed at runtime (for example `angles[idx]`, with `idx` taken from a measurement result).\n", "\n", "Both are supported, but dynamic indexing is more general and usually a bit more demanding for lowering/compilation. The examples below illustrate both patterns." ] }, { "cell_type": "code", "execution_count": 8, "id": "27ac7ebd", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Static indexing (single run): False\n", "Static indexing (10 shots): [False, False, False, False, False, False, False, False, False, False]\n" ] } ], "source": [ "import numpy as np\n", "import jax.numpy as jnp\n", "from qrisp_cudaq import FixedShapeNDArray\n", "\n", "# Static indexing example: access angles[0], where the index is known at compile-time.\n", "@cudaq_kernel\n", "def static_index_kernel(angles: FixedShapeNDArray[float, 3]):\n", " qv = QuantumFloat(1)\n", " rx(angles[0], qv[0])\n", " return measure(qv[0])\n", "\n", "angles_static = np.array([0.0, 1.57, 3.14])\n", "result_static = static_index_kernel(angles_static)\n", "print(\"Static indexing (single run):\", result_static)\n", "\n", "shots_static = cudaq.run(static_index_kernel, angles_static, shots_count=10)\n", "print(\"Static indexing (10 shots):\", shots_static)" ] }, { "cell_type": "markdown", "id": "742b9a0e", "metadata": {}, "source": [ "In the previous cell, `angles[0]` is fixed at compile time. Initially, the qubit is in $\\ket{0}$ and `rx(0.0)` leaves it there, so the result is `False`. \n", "\n", "Now let's switch to **dynamic indexing**: we compute the index at runtime from a measurement result and then read `angles[idx]`. The Hadamard puts the qubit into an equal superposition, so the measurement gives `idx = 0` or `idx = 1` with equal probability, and the kernel applies either `angles[0]` or `angles[1]`." ] }, { "cell_type": "code", "execution_count": 9, "id": "24896996", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Dynamic indexing (single run): True\n", "Dynamic indexing (10 shots): [False, False, False, False, False, False, False, False, False, False]\n" ] } ], "source": [ "# Dynamic indexing example: index is computed at runtime.\n", "@cudaq_kernel\n", "def dynamic_index_kernel(angles: FixedShapeNDArray[float, 3]):\n", " qv = QuantumFloat(1)\n", " h(qv)\n", " idx = jnp.int32(measure(qv[0]))\n", " rx(angles[idx], qv[0])\n", " return measure(qv[0])\n", "\n", "angles_dynamic = np.array([0.0, 1.57, 3.14])\n", "result_dynamic = dynamic_index_kernel(angles_dynamic)\n", "print(\"Dynamic indexing (single run):\", result_dynamic)\n", "\n", "shots_dynamic = cudaq.run(dynamic_index_kernel, angles_dynamic, shots_count=10)\n", "print(\"Dynamic indexing (10 shots):\", shots_dynamic)" ] }, { "cell_type": "markdown", "id": "b1cd45c7", "metadata": {}, "source": [ "If `idx` is 0, the qubit is in $\\ket{0}$ and `rx(0.0)` leaves it there, so the result is `False`. If `idx` is 1, the qubit is in $\\ket{1}$ and `rx(1.57)` (about $\\pi/2$) puts it into an equal superposition, so the result is `True` half of the time. Overall, about a quarter of the shots return `True`." ] }, { "cell_type": "markdown", "id": "21da61ce", "metadata": {}, "source": [ "### Static vs Dynamic Arrays in CUDA-Q Kernels\n", "\n", "This topic is important because array behavior changes depending on where the array comes from.\n", "\n", "- **Static array**: created inside the kernel body (for example via `np.array([...])`).\n", "- **Dynamic array parameter**: passed into the kernel and annotated as `FixedShapeNDArray[dtype, n]` or created inside the kernel (i.e., `jnp.array([...])`).\n", "\n", "Practical rule of thumb:\n", "\n", "- Static arrays defined inside the kernel support standard array transformations (for example `np.sin`).\n", "- Dynamic array parameters are primarily meant for indexed access (`angles[i]`) and passing values into quantum operations. Transformations on dynamic arrays inside a kernel (for example `jnp.sin`) are not supported at this point.\n", "\n", "The three examples below demonstrate these patterns in executable form." ] }, { "cell_type": "code", "execution_count": 10, "id": "1fbfd378", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Static array in-kernel (10 shots): [0.0, 0.0, 0.0, 3.0, 3.0, 0.0, 1.0, 0.0, 1.0, 0.0]\n" ] } ], "source": [ "# Example 1: static array created inside the kernel\n", "@cudaq_kernel\n", "def static_array_kernel():\n", " qv = QuantumFloat(5)\n", "\n", " angles = np.array([1.57, 0.78, 0.39, 0.25, 0.12])\n", " angles = np.sin(angles)\n", "\n", " # Static array with static loop.\n", " for i in range(5):\n", " ry(angles[i], qv[i])\n", "\n", " return measure(qv)\n", "\n", "print(\"Static array in-kernel (10 shots):\", cudaq.run(static_array_kernel, shots_count=10))" ] }, { "cell_type": "code", "execution_count": 11, "id": "b8afc961", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Dynamic array parameter (10 shots): [0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 3.0, 0.0, 1.0, 1.0]\n" ] } ], "source": [ "# Example 2: dynamic array passed as kernel parameter\n", "@cudaq_kernel\n", "def dynamic_array_kernel(angles: FixedShapeNDArray[float, 5]):\n", " qv = QuantumFloat(5)\n", "\n", " # Dynamic parameter array accessed element-wise.\n", " for i in jrange(5):\n", " ry(angles[i], qv[i])\n", "\n", " return measure(qv)\n", "\n", "angles = np.array([1.57, 0.78, 0.39, 0.25, 0.12])\n", "print(\"Dynamic array parameter (10 shots):\", cudaq.run(dynamic_array_kernel, angles, shots_count=10))" ] }, { "cell_type": "code", "execution_count": 12, "id": "9f74f1bd", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Error defining kernel: Failed to compile Qrisp function 'dynamic_array_bad_kernel' to CUDA-Q: This @cudaq_kernel function performs arithmetic on a dynamic (traced) array, which CUDA-Q cannot compile.\n", "\n", "Static Python arrays created inside the kernel are supported. JAX arrays are dynamic during tracing, and every array passed as a kernel argument is dynamic. CUDA-Q does not support arithmetic on dynamic arrays inside a kernel.\n", "\n", "To fix this:\n", "- Move JAX array calculations outside the @cudaq_kernel function and pass scalar results into it.\n", "- When passing an array to the kernel, access individual elements instead of combining or adding the whole array.\n", "- Use scalar values or static Python arrays for calculations that stay inside the kernel.\n", "\n", "\n" ] } ], "source": [ "# Example 3: unsupported dynamic array operations inside the kernel (should fail to compile)\n", "try:\n", " @cudaq_kernel\n", " def dynamic_array_bad_kernel():\n", " qv = QuantumFloat(5)\n", "\n", " # Operations like jnp.sin on dynamic arrays are not currently supported inside the kernel.\n", " angles = jnp.array([1.57, 0.78, 0.39, 0.25, 0.12])\n", " angles = jnp.sin(angles)\n", "\n", " # Dynamic array with static loop.\n", " for i in range(5):\n", " ry(angles[i], qv[i])\n", "\n", " return measure(qv)\n", "except RuntimeError as e:\n", " print(\"Error defining kernel:\", e)" ] }, { "cell_type": "markdown", "id": "84e1bdb8", "metadata": {}, "source": [ "## Variational Algorithms: QAOA\n", "\n", "So far, classical logic ran *inside* the kernel, during each shot. Variational algorithms are hybrid in a different way: the kernel makes no decisions based on measurement results. Instead, a classical optimizer running in ordinary Python calls the kernel again and again with new parameters and steers them toward better results. The hybrid loop runs *between* kernel executions, and the kernel parameters from the previous section are what make it possible: Qrisp compiles the kernel once, and each call only passes in new parameter values.\n", "\n", "The best-known example is the **Quantum Approximate Optimization Algorithm (QAOA)**. We use it for **MaxCut**: split the nodes of a graph into two groups so that as many edges as possible connect nodes in different groups. Each qubit stands for one node, and its measured value says which group the node belongs to. QAOA starts in a superposition of all possible splits and alternates two operations for a fixed number of layers:\n", "\n", "- the **cost operator** gives each split a phase that depends on how many edges it cuts, scaled by an angle $\\gamma$;\n", "- the **mixer** rotates every qubit by an angle $\\beta$, moving amplitude between different splits.\n", "\n", "The angles $\\gamma$ and $\\beta$ of each layer are the kernel's parameters, passed in as a `FixedShapeNDArray`. A SciPy optimizer adjusts them to maximize the average cut of the measured bitstrings.\n", "Since the kernel makes no decisions based on measurement results, it can use `execution_mode=\"sample\"`: `cudaq.sample` returns a histogram of the measured bitstrings. This is the execution mode most hardware backends support.\n" ] }, { "cell_type": "code", "execution_count": 13, "id": "094beeb5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Starting optimization... Initial expected cut: 5.370\n", "\n", "Optimization finished! ✨\n", "Final expected cut: 7.300\n", "Best bitstring found: 1010101\n", "Max cut value achieved: 9\n", "Optimized parameters: [1.552 1.304 2.182 1.99 0.346 0.331]\n" ] } ], "source": [ "from scipy.optimize import minimize\n", "import numpy as np\n", "\n", "# 1. Define the MaxCut Instance\n", "EDGES = [(0, 1), (1, 2), (2, 3), (3, 0), (0, 2), (1, 3), (1, 4), (2, 5), (3, 6), (4, 5), (5, 6), (6, 4)]\n", "N = 7\n", "\n", "# 2. Define the QAOA operators and the number of layers\n", "num_layers = 3\n", "\n", "def apply_cost_operator(qv, gamma):\n", " for u, v in EDGES:\n", " cx(qv[u], qv[v])\n", " rz(2 * gamma, qv[v])\n", " cx(qv[u], qv[v])\n", "\n", "def apply_mixer(qv, beta):\n", " for i in range(N):\n", " rx(2 * beta, qv[i])\n", "\n", "# 3. Define the Qrisp CUDA-Q Kernel (sample mode: no return value)\n", "@cudaq_kernel(execution_mode=\"sample\")\n", "def variational_maxcut_kernel(params: FixedShapeNDArray[float, 2 * num_layers]):\n", " qv = QuantumVariable(N)\n", "\n", " # Start in the uniform superposition\n", " h(qv)\n", "\n", " # Apply multiple QAOA layers\n", " for i in jrange(num_layers):\n", " apply_cost_operator(qv, params[i]) # Gammas\n", " apply_mixer(qv, params[i + num_layers]) # Betas\n", "\n", " measure(qv)\n", "\n", "# 4. Classical Post-Processing\n", "def calculate_cut_value(bitstring: str) -> int:\n", " \"\"\"Calculates the cut size of a measured bitstring (character i is node i).\"\"\"\n", " return sum(bitstring[u] != bitstring[v] for u, v in EDGES)\n", "\n", "def average_cut(counts) -> float:\n", " \"\"\"Average cut value over all shots of a sampled histogram.\"\"\"\n", " shots = sum(count for _, count in counts.items())\n", " return sum(calculate_cut_value(b) * count for b, count in counts.items()) / shots\n", "\n", "def objective_function(params: np.ndarray) -> float:\n", " \"\"\"Samples the quantum kernel and evaluates the expected (average) cut.\"\"\"\n", " # Sample the circuit: a histogram of bitstrings and their counts\n", " counts = cudaq.sample(variational_maxcut_kernel, params, shots_count=100)\n", "\n", " # SciPy minimizes, so we return the negative expected cut\n", " return -average_cut(counts)\n", "\n", "# 5. Hybrid Optimization Loop\n", "np.random.seed(42)\n", "initial_params = np.random.uniform(0.0, np.pi/2, size=2 * num_layers)\n", "\n", "print(f\"Starting optimization... Initial expected cut: {-objective_function(initial_params):.3f}\")\n", "\n", "opt_result = minimize(\n", " objective_function,\n", " initial_params,\n", " method=\"COBYLA\",\n", " options={\"maxiter\": 40}\n", ")\n", "\n", "# 6. Evaluate the Best Result\n", "best_params = opt_result.x\n", "final_counts = cudaq.sample(variational_maxcut_kernel, best_params, shots_count=100)\n", "\n", "# Find the best bitstring in our final sampled distribution\n", "best_bitstring = max((b for b, _ in final_counts.items()), key=calculate_cut_value)\n", "\n", "print(\"\\nOptimization finished! ✨\")\n", "print(f\"Final expected cut: {average_cut(final_counts):.3f}\")\n", "print(f\"Best bitstring found: {best_bitstring}\")\n", "print(f\"Max cut value achieved: {calculate_cut_value(best_bitstring)}\")\n", "print(f\"Optimized parameters: {np.round(best_params, 3)}\")" ] }, { "cell_type": "markdown", "id": "07bdb0e4", "metadata": {}, "source": [ "During the optimization, the average cut of the measured bitstrings improves. It starts close to 6, which is what random guessing achieves on this graph: a random split cuts half of the 12 edges on average. The best sampled bitstring describes a split of the nodes into two groups, where the $i$-th character from the left stands for node $i$. A cut of 9 is the best possible for this graph.\n", "\n", "Keep in mind that this instance is deliberately tiny. With 7 nodes there are only $2^7 = 128$ possible splits, so even random guessing with 100 shots almost certainly finds an optimal one, and trying all splits classically is instant. QAOA becomes interesting for larger graphs, where the number of splits, $2^N$, far exceeds any number of shots, and the algorithm's job is to concentrate the measurement probability on good cuts. The workflow shown here stays the same; only the graph changes.\n" ] }, { "cell_type": "markdown", "id": "825011f8", "metadata": {}, "source": [ "## Recap\n", "\n", "So far we've covered a lot of ground. Let's gather it into a list:\n", "\n", "- **Bell state & basic kernels**: defining `@cudaq_kernel` functions and running them with `cudaq.run`, or with `cudaq.sample` for kernels compiled with `execution_mode=\"sample\"`.\n", "- **Real-time hybrid logic**: classical logic that runs inside the kernel during every shot, such as feed-forward corrections in quantum teleportation, a `q_while_loop` whose length depends on measurement results, and iterative phase estimation, where each round of a `q_fori_loop` applies a phase correction computed from the earlier measurements.\n", "- **Kernel basics**: multiple return values, scalar (`int`, `float`, `bool`) and array (`FixedShapeNDArray`) parameters, static vs dynamic indexing, and which array operations CUDA-Q can compile.\n", "- **Variational hybrid loops**: QAOA for MaxCut, where a SciPy optimizer calls a parametrized kernel from ordinary Python.\n", "\n", "All of this runs natively on CUDA-Q. But the real question is: can we cover even more ground?\n", "\n", "[Block encodings](https://www.qrisp.eu/general/tutorial/BE_tutorial/BE_vol1.html) are one of the most powerful and mathematically sophisticated tools in Qrisp. They represent linear operators as the top-left block of a unitary, enabling polynomial transformations of spectra, matrix inversion, Hamiltonian simulation, and much more. If the CUDA-Q interface can handle block encodings, it can handle essentially anything Qrisp throws at it.\n", "\n", "So, can this interface run block encodings? Let's find out!\n" ] }, { "cell_type": "markdown", "id": "9d807034", "metadata": {}, "source": [ "## Final Example: Block Encodings\n", "\n", "A [block encoding](https://arxiv.org/abs/1806.01838) of a matrix $A$ is a unitary $U$ such that\n", "\n", "$$ A = (\\langle 0 | \\otimes I) U (|0 \\rangle \\otimes I) $$\n", "\n", "i.e., $A$ appears as the top-left block of $U$ when the ancilla variables are projected onto $|0\\rangle$. This construction is the foundation of quantum singular value transformation (QSVT) and enables a broad class of quantum algorithms: Hamiltonian simulation, quantum linear systems solvers, spectral filtering, and more.\n", "\n", "Qrisp's [BlockEncoding](../../reference/Block%20Encodings/BlockEncoding.rst) class lets you construct block encodings directly from operator expressions and compose them with `+`, `-`, `@`, `.poly()`, `.inv()`, `.sim()`, and so on, all **outside** the kernel, since those compositions involve classical arithmetic on traced JAX arrays, which is not supported inside the tracing context.\n", "\n", "Inside the kernel, you call `BE.apply(operand)` to apply the block encoding and obtain the ancilla variables. Because the encoding is only exact up to post-selection on the ancillas being in $|0\\rangle$, the kernel returns **both** the operand measurement and the ancilla measurements. The classical post-processing then filters for the `ancilla == 0` outcomes." ] }, { "cell_type": "code", "execution_count": 14, "id": "b1e2225f", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Post-selected distribution: {0: 0.8989898989898989, 8: 0.0404040404040404, 2: 0.0404040404040404, 1: 0.0202020202020202}\n" ] } ], "source": [ "from collections import Counter\n", "from qrisp.block_encodings import BlockEncoding\n", "from qrisp.operators import X, Y, Z\n", "\n", "\n", "def post_selection(res_dict):\n", " # Post-selection on ancillas being in the |0> state\n", " filtered_dict = {k[0]: p for k, p in res_dict.items() \\\n", " if all(x == 0 for x in k[1:])}\n", " success_prob = sum(filtered_dict.values())\n", " filtered_dict = {k: p / success_prob for k, p in filtered_dict.items()}\n", " return filtered_dict\n", "\n", "\n", "# Define an Ising Hamiltonian H with a transverse field and its block encoding,\n", "# then use the block encoding for Hamiltonian simulation: BE.sim(t, N) block-encodes\n", "# the time evolution e^{-itH}. All of this happens outside the kernel.\n", "H = (1/7)*(Z(0)*Z(1) + Z(1)*Z(2) + Z(2)*Z(3) + X(0) + X(1) + X(2) + X(3))\n", "BE = BlockEncoding.from_operator(H)\n", "BE_sim = BE.sim(t=1.0, N=3)\n", "\n", "@cudaq_kernel\n", "def main():\n", "\n", " operand = QuantumVariable(4)\n", " ancs = BE_sim.apply(operand)\n", " return measure(operand), *[measure(anc) for anc in ancs]\n", "\n", "results = cudaq.run(main, shots_count=100)\n", "counts = Counter(results)\n", "res_dict = {outcome: count / len(results) for outcome, count in counts.items()}\n", "filtered_dict = post_selection(res_dict)\n", "print(\"Post-selected distribution:\", filtered_dict)" ] }, { "cell_type": "markdown", "id": "e10d71cd", "metadata": {}, "source": [ "This example shows how to apply the time evolution $e^{-iH}$ for the Ising Hamiltonian $H$ to a variable in the initial state $\\ket{\\psi_0}=\\ket{0}$. The block encoding is applied successfully if the ancilla variables are measured in state $\\ket{0}$. Since $e^{-iH}$ itself is unitary, the success probability will be close to 1, up to the approximation error. The post-selected distribution shows the measurement distribution of the state $\\ket{\\psi}=e^{-iH}\\ket{\\psi_0}$. After the short evolution time, the state is still mostly $\\ket{0}$, but the transverse field has started to flip individual spins, so the outcomes 1, 2, 4 and 8 (one flipped spin each) appear in a few percent of the shots." ] }, { "cell_type": "markdown", "id": "ce962c01", "metadata": {}, "source": [ "## Epilogue\n", "\n", "You made it! By now, you have seen how the Qrisp CUDA-Q interface bridges the gap between high-level quantum abstractions and hardware-accelerated execution. We started with a humble Bell state, teleported a qubit, estimated the energy of a quantum magnet, navigated the nuances of static versus dynamic arrays, built a variational QAOA loop, and even wrangled block encodings.\n", "\n", "All these features combine into flexible quantum-classical workflows, and the true power of this integration lies in what you build next. Whether you are scaling up variational pipelines, experimenting with dynamic-circuit feed-forward, or diving deep into quantum singular value transformation, the heavy lifting is handled.\n", "\n", "Happy coding, and may your quantum states remain coherent!\n" ] }, { "cell_type": "markdown", "id": "73cf3073", "metadata": {}, "source": [ "---\n", "\n", "## Tips, Gotchas, and Debugging\n", "\n", "- **Kernel parameters**: Use `FixedShapeNDArray` for NumPy arrays, and regular Python types for scalars.\n", "- **Hybrid control**: You can use `measure()`, `control()`, `q_cond()`, `q_while_loop()`, and `q_fori_loop()` to build hybrid quantum-classical logic.\n", "- **Debugging**: If you get errors, check your kernel signature and make sure all types are supported by CUDA-Q.\n", "- **Limitations**: Some advanced Qrisp features may not yet be supported inside kernels. For that reason, keep kernels simple and explicit for best results.\n", "\n", "If you run into trouble, open an issue or ask other Qrisp devs for help and/or insights on the [Eclipse Qrisp Discord server](https://discord.gg/v5np7DeBaq).\n", "\n", "---" ] } ], "metadata": { "kernelspec": { "display_name": "qrisp_dev", "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.11.15" } }, "nbformat": 4, "nbformat_minor": 5 }