Debugging with `state_result` statements¶

Download this notebook - state_results.ipynb

In this example we will demonstrate the use of state_result statements for surfacing the simulator state while running a Guppy program with Selene.

from guppylang import guppy
from guppylang.std.builtins import array
from guppylang.std.debug import state_result
from guppylang.std.quantum import qubit, discard, discard_array, cx, h, measure, x

from hugr.qsystem.result import QsysResult

from selene_sim import build, Quest

import numpy as np

np.set_printoptions(precision=4, suppress=True, linewidth=120)

Consider the following Guppy program where we create a 3-qubit GHZ state in two different ways that should be unitarily equivalent.

@guppy
def ghz_state(q: array[qubit, 3]) -> None:
    h(q[0])
    cx(q[0], q[1])
    # Apply CX to second and third qubit.
    cx(q[1], q[2])


@guppy
def ghz_state_alternative(q: array[qubit, 3]) -> None:
    h(q[0])
    cx(q[0], q[1])
    # Apply CX to first and third qubit.
    cx(q[0], q[2])

In order to check whether the resulting state vectors are equal, we can use the state_result function. As inputs it takes a tag (just as result statements do) and then either one array of qubits or a variable number of single qubit arguments. The order of arguments later determines the order of state vector outputs.

@guppy
def main() -> None:
    qs = array(qubit() for _ in range(3))
    ghz_state(qs)
    state_result("ghz", qs)
    discard_array(qs)

    qs = array(qubit() for _ in range(3))
    ghz_state_alternative(qs)
    state_result("ghz_alternative", qs)
    discard_array(qs)

Let’s now extract the requested state from the simulator. This can be done by calling extract_states_dict on the simulator plugin for each shot to get a dictionary mapping tags to states, and then retrieving the results from the state using get_single_state.

Here our ghz_state and ghz_state_alternative functions deterministically prepare a single quantum state.

Note that this is currently only possible using the Quest simulator plugin with the default SimpleRuntime and IdealErrorModel.

shots = main.emulator(n_qubits=3).statevector_sim().with_seed(42).run()

for states in shots.partial_state_dicts():
    state_vector1 = states["ghz"].as_single_state()
    print(f"State `ghz1`: {state_vector1}")

    state_vector2 = states["ghz_alternative"].as_single_state()
    print(f"State `ghz2`: {state_vector2}")

    assert np.allclose(state_vector1, state_vector2)

State `ghz1`: [0.7071+0.j 0.    +0.j 0.    +0.j 0.    +0.j 0.    +0.j 0.    +0.j 0.    +0.j 0.7071+0.j]

State `ghz2`: [0.7071+0.j 0.    +0.j 0.    +0.j 0.    +0.j 0.    +0.j 0.    +0.j 0.    +0.j 0.7071+0.j]

We can see that the state vector is the same for both GHZ functions!

It is important to bear in mind that the state result is simply the internal state of the simulator plugin, so the output at any given point in the middle of the program could be very different on hardware, where various optimisations and different scheduling likely apply.

If you chose to request a subset of the qubits making up the overall state, you should also be aware that qubits in this subset could be entangled with other qubits in the system, so the qubits not requested in the state result are traced out to provide a probabilistic set of state vectors. For example, if we only request the state of any two out of the three qubits in the GHZ state, we expect to see either \(\ket{00}\) or \(\ket{11}\) with a probability of \(\frac{1}{2}\).

@guppy
def main() -> None:
    qs = array(qubit() for _ in range(3))
    ghz_state(qs)
    state_result("subset01", qs[0], qs[1])
    state_result("subset02", qs[0], qs[2])
    state_result("subset12", qs[1], qs[2])
    discard_array(qs)


runner = build(main.compile())

shots = main.emulator(n_qubits=4).statevector_sim().with_seed(42).run()


for states in shots.partial_state_dicts():
    state01 = states["subset01"].state_distribution()
    state02 = states["subset02"].state_distribution()
    state12 = states["subset12"].state_distribution()

    print(
        f"{state01[0].state} -> {round(state01[0].probability, 2)}, {state01[1].state} -> {round(state01[1].probability, 2)}"
    )
    print(
        f"{state02[0].state} -> {round(state02[0].probability, 2)}, {state02[1].state} -> {round(state02[1].probability, 2)}"
    )
    print(
        f"{state12[0].state} -> {round(state12[0].probability, 2)}, {state12[1].state} -> {round(state12[1].probability, 2)}"
    )

[1.+0.j 0.+0.j 0.+0.j 0.+0.j] -> 0.5, [0.+0.j 0.+0.j 0.+0.j 1.+0.j] -> 0.5
[1.+0.j 0.+0.j 0.+0.j 0.+0.j] -> 0.5, [0.+0.j 0.+0.j 0.+0.j 1.+0.j] -> 0.5
[1.+0.j 0.+0.j 0.+0.j 0.+0.j] -> 0.5, [0.+0.j 0.+0.j 0.+0.j 1.+0.j] -> 0.5

Also keep in mind that state results can vary between shots depending on control flow.

@guppy
def main() -> None:
    q = qubit()
    h(q)
    q1 = qubit()
    q2 = qubit()
    x(q2)
    if measure(q):
        state_result("state", q2)
    else:
        state_result("state", q1)
    discard(q1)
    discard(q2)

shots = main.emulator(n_qubits=4).statevector_sim().with_seed(2).with_shots(100).run()


state_counts = {0: 0, 1: 0}
for state_dict in shots.partial_state_dicts():
    state = state_dict["state"].state_distribution()
    s = state[0].state
    if np.allclose(s, [0, 1]):
        state_counts[1] += 1
    elif np.allclose(s, [1, 0]):
        state_counts[0] += 1

print(state_counts)

{0: 58, 1: 42}

Debugging with state_result statements¶

Debugging with `state_result` statements¶