Add state_result support for Stim

[jake-arkinstall]

state_result is currently only supported with the QuEST plugin.

This is the logic of state printing:

• A user can call state_result("some_label", qubit, another_qubit, some_qubit_list)
• In this case we call the 'qubit, another_qubit, some_qubit_list' list the "requested qubits".
• This is generally a subset of the available qubits.
• They need to be able to extract meaningful state information from this

The caveat: the qubits they request are generally not the full qubit register, so the 'state' of the qubits they request is generally a mixed state.

We could do with the equivalent functionality for the Stim plugin so that users can inspect the state of larger registers, either with a small number of requested qubits (in which case a statevector/density matrix representations are reasonable) or with a large number of requested qubits (in which case equivalent tableau representations are reasonable)

[jake-arkinstall]

Here is how it is handled in the QuEST plugin:

• The state of the full register is written to a file in binary format, including the number of qubits, the number of requested qubits, and the indices of requested qubits
• The result stream is provided with the path to that file
• The user can convert this into a SeleneQuestState, which opens the file and parses the data
• The state member stores the full state, and an array of indices of requested qubits is stored too
• This allows users to access the full state (not typically meaningful) but also a partial density matrix and a list of states with their associated probabilities.
• That's done by tensor manipulations in numpy

To get similar functionality in the Stim plugin:

• We should dump the stabiliser tableau to the binary file, along with metadata like the requested qubits
• Should this be read directly into the Stim python package? Should we add that as a dependency?
• Alternatively we could perform tableau manipulations directly in order to confine the state to requested qubits:
  • I believe that it is feasible to get a density matrix representation if the number of requested qubits is small, but for a large number of qubits this require unrealistic amounts of processing and storage.
  • Is there a mixed state representation of tableaus?
  • A stabiliser-native alternative to producing a list of states with probabilities is to produce a list of tableaus with their probabilities. This should be viable, but what's the optimal approach? A naive approach is to postselect on every combination of unrequested qubit outcomes, but that would be enormously inefficient.

[PabloAndresCQ]

Is there a mixed state representation of tableaus?

Yes. You can get an efficient representation of "a list of tableaus". Each of the tableaus would have the same m Pauli strings of m qubits that the user selected, the different tableaus just differ on whether the state they represent is a +1 or -1 eigenvector of the each of the individual m Pauli strings. You can represent this by including a boolean function per Pauli string which indicates if it's +1 or -1. The boolean function has n-m variables (where n is the number of qubits in the "full state") and it's just made up of XOR operations. You then know all 2^{n-m} stabiliser states that would come from postselecting the discarded qubits by choosing any combination of values for those boolean variables (in fact, there's only at most 2^m distinct states).

Calculating these boolean functions is straightforward, you just need to do some partial Gaussian elimination on the original n qubit tableau to put n-m of the rows in a special form, keep track of the row-sums you applied to achieve this (which gives you the boolean functions on each of the m remaining rows) and finally discard the n-m rows. I can explain it in more detail at some point.

I expect this representation will be easy to interpret by users that are fluent on the stabiliser formalism. I also think that if m is small (and no matter how large n is), this can be efficiently converted to the list-of-statevectors-with-probability representation that we already had for QuEST: there can only be 2^m different tableaus no matter how large n is, and figuring out whether two different assignments of the n variables give the same tableau should just be polynomial on n (it's just a matter of solving linear systems of equations).