Sample from wavefunction has wrong bit order

[jj-curious]

Hi,

I've noticed that the bit order is wrong when sampling from a wavefunction with the sample_bitstrings() method. It returns bitstrings in little endian order, which is different from using a QPU / QVM (big endian order).

How to Reproduce

Code Snippet

import pyquil

p = pyquil.Program()
p += pyquil.gates.X(0)
p += pyquil.gates.X(1)
p += pyquil.gates.I(2)

sim = pyquil.api.WavefunctionSimulator()
wf = sim.wavefunction(p)
bitstrings = wf.sample_bitstrings(5)
print(bitstrings)

Error Output

This should actually return [1 1 0] (using a QVM does indeed return [1 1 0]).

[[0 1 1]
 [0 1 1]
 [0 1 1]
 [0 1 1]
 [0 1 1]]

Proposed solution

This can be fixed by modifying the sample_bitstrings() method from https://github.com/rigetti/pyquil/blob/master/pyquil/wavefunction.py in line 200-210.

For example. simply reversing the bitstrings in the return statement using [:, ::-1].

    def sample_bitstrings(self, n_samples: int) -> np.ndarray:
        """
        Sample bitstrings from the distribution defined by the wavefunction.
        :param n_samples: The number of bitstrings to sample
        :return: An array of shape (n_samples, n_qubits)
        """
        possible_bitstrings = np.array(list(itertools.product((0, 1), repeat=len(self))))
        inds = np.random.choice(2 ** len(self), n_samples, p=self.probabilities())
        bitstrings = possible_bitstrings[inds, :]
        return bitstrings[:, ::-1]   # Here: need to reverse the bitstrings

[notmgsk]

Yeah, that's unfortunate. Unfortunate because it's not really something we can go back and change now -- anybody using that code will expect the current ordering. After a change, that would break. :(

What we can do is make it clear in the documentation. You're welcome to make a pull request :)

[mhodson-rigetti]

For completeness, the following additional code gives the (different) QVM result:

qvm = pyquil.get_qc("3q-qvm")
ro = p.declare("ro", "BIT", 3)
p += pyquil.gates.MEASURE(0, ro[0])
p += pyquil.gates.MEASURE(1, ro[1])
p += pyquil.gates.MEASURE(2, ro[2])
p.wrap_in_numshots_loop(5)
bitstrings = qvm.run(p)
print(bitstrings)

[mhodson-rigetti]

wavefunction_fast_sample_bitstrings.zip

This notebook provides an alternative solution that (a) dramatically improves the speed of generating samples (more than a factor of 10 in some cases); (b) provides an "as_qam" flag that corrects the bit ordering issue raised here, but optionally (not breaking backward compatibility); and (c) also controls the random seed, and so addresses https://github.com/rigetti/pyquil/issues/1272.

The function in this notebook operates on the Wavefunction object, but could easily be refactored to replace the intrinsic method. Recommend that be done and tests added around it to further verify and check for / prevent future regression.