Measuring multiple qubit arrays

[perlinm]

I need to save the state vector of qubits that organized into several arrays, many_arrays: array[array[qubit, N], M]. The number of arrays and their sizes are both known at compile time, but may otherwise vary. I was hoping to do something like

state_result("state_vector", array(qubit for one_array in many_arrays for qubit in one_array))

but this throws the error

Nested array comprehensions are not supported

How can I save the state vector of all qubits?

[GBoydQ]

Hi, I have worked around similar issues to this in the past before this is supported, using some knowledge of the internals of the SeleneQuestState object that is produced from state_results although what might be easiest will depend on your specific scenario.

Essentially SeleneQuestStates always store the full statevector of the simulation, and the argument in state_result just determines a list of specified qubits that are then used to get partial state information if you ask for the density matrix (or pure state assuming no entanglement), so you can manipulate this specified_qubits member after your simulation to get your desired results.

1. The easiest scenario is that these qubits in many_arrays are all the qubits that have been used in your simulation (including qubits that were discarded earlier in the program). In this case, any state result will already contain the information you want in SeleneQuestState.state (note of warning, the QuEST endianness convention of the internal simulation is opposite to guppy's usual convention). (You could also set specified_qubits = list(range(N*M) see below)

2. If this is not the case, you can create multiple state results and manually join their specified qubits e.g.

@guppy
def main()->None:
    reg1 = array(qubit() for _ in range(3))
    reg2 = array(qubit() for _ in range(2))
    unimportant_qubit = qubit()
    state_result("reg1", reg1)
    state_result("reg2", reg2)

res = main.emulator(6).run()
states = Quest.extract_states_dict(res.results[0].entries)

reg1_bits = states["reg1"].specified_qubits
reg2_bits = states["reg2"].specified_qubits

both_reg_state = states["index"]
both_reg_state.specified_qubits = reg1_bits + reg2_bits

this might be slightly more difficult in your case as you may not be able to use string literals for the state_results as you loop over M.

[GBoydQ]

My bad, I didn't read your code properly

states = Quest.extract_states_dict(result.results[0].entries)
state_vector = states["state_vector_0"].state

[GBoydQ]

where Quest is
from selene_sim import Quest

[perlinm]

Great, thank you! This works. It is still a hack, though, so I'll leave the issue open.

[perlinm]

Ah, I have one more issue: I do not now which qubit array lives where in memory, because mem.mem_swap is used in the program. So I need to do something closer to the original example here in terms of combining specified_qubits, but now I am back to figuring out how to do this for N registers. Unfortunately, this throws a Guppy error:

for ii in range(comptime(num_qubit_arrays)):
    state_result(f"state_vector_{ii}", many_arrays[ii])

because I think Guppy is unhappy with not knowing f"state_vector_{ii}" at compile time.

Also, it's not clear but it seems like multiple state_result commands may store redundant copies of the state vector, which would not be ideal.

[perlinm]

Okay, so my current hack looks something like this. The Guppy program has

        if comptime(num_qubit_arrays) > 0:
            state_result("state_vector_0", qubit_arrays[0])
        if comptime(num_qubit_arrays) > 1:
            state_result("state_vector_1", qubit_arrays[1])
        if comptime(num_qubit_arrays) > 2:
            state_result("state_vector_2", qubit_arrays[2])
        if comptime(num_qubit_arrays) > 3:
            state_result("state_vector_3", qubit_arrays[3])
        if comptime(num_qubit_arrays) > 4:
            state_result("state_vector_4", qubit_arrays[4])
        if comptime(num_qubit_arrays) > 5:
            panic("Saving the state vector of more than 5 qubit arrays is unsupported")

In theory, this should work to extract the state vector:

        sim_result = program.emulator(n_qubits=num_physical_qubits).statevector_sim().run()
        states = selene_sim.Quest.extract_states_dict(sim_result.results[0].entries)
        state = states["state_vector_0"]
        for ii in range(1, num_qubit_arrays):
            state.specified_qubits += states[f"state_vector_{ii}"].specified_qubits
        state_vector = state.get_single_state()

But in practice I find that if state.specified_qubits is too long, the last step of state.get_single_state() is either extremely slow or freezes altogether (I am guessing due to the slowness of checking that the requested qubits are in a pure state), so instead I do the following:

        num_qubits = ...
        num_qubit_arrays = ...
        sim_result = program.emulator(n_qubits=num_qubits).statevector_sim().run()
        states = selene_sim.Quest.extract_states_dict(sim_result.results[0].entries)
        data_qubits = [
            data_qubit
            for qubit_array_index in range(num_qubit_arrays)
            for data_qubit in states[f"state_vector_{qubit_array_index}"].specified_qubits
        ]
        tensor = states["state_vector_0"].state.reshape([2] * num_qubits).T
        tensor = np.moveaxis(tensor, data_qubits, range(len(data_qubits)))
        tensor = tensor.reshape([2 ** len(data_qubits), -1])
        for state_vector in tensor.T:
            if not np.allclose(state_vector, 0):
                state_vector = state_vector / np.linalg.norm(state_vector)
                break

Of course, this assumes that the joint state of the qubit arrays is a pure state.