`run_and_measure` with `define_noisy_gate` gives biased results

Using define_noisy_gate to simulate a channel, I get expected results if I use run but not when I use run_and_measure. More specifically, if I run

%matplotlib inline
import matplotlib.pyplot as plt
from pyquil import get_qc, Program
import numpy as np


def kraus_ops_bit_flip(prob):
    # define flip (X) and not flip (I) Kraus operators
    I_ = np.sqrt(1 - prob) * np.array([[1, 0], [0, 1]])
    X_ = np.sqrt(prob) * np.array([[0, 1], [1, 0]])
    return [I_, X_]

def random_unitary(n):
    # draw complex matrix from Ginibre ensemble
    z = np.random.randn(n, n) + 1j * np.random.randn(n, n)
    # QR decompose this complex matrix
    q, r = np.linalg.qr(z)
    # make this decomposition unique
    d = np.diagonal(r)
    l = np.diag(d) / np.abs(d)
    return np.matmul(q, l)

# pick probability
prob = 0.2
# noisy program
p = Program()
p.defgate("DummyGate", random_unitary(2))
p += ("DummyGate", 0)
p.define_noisy_gate("DummyGate", [0], kraus_ops_bit_flip(prob))

qc = get_qc('1q-qvm')

num_expts = 1000
num_shots = 1000
results = np.zeros(num_expts)

for i in range(num_expts):
    results_tmp = qc.run_and_measure(p, trials=num_shots)
    results[i] = np.count_nonzero(results_tmp[0]) / num_shots

    
plt.figure(figsize=(10, 8))
plt.hist(results, bins=50)
plt.show()

I get the following plot, which should really be centered at 0.2 but is instead centered on some other value:

If instead I run the following code using run:

%matplotlib inline
import matplotlib.pyplot as plt
from pyquil import get_qc, Program
from pyquil.gates import MEASURE
import numpy as np


def kraus_ops_bit_flip(prob):
    # define flip (X) and not flip (I) Kraus operators
    I_ = np.sqrt(1 - prob) * np.array([[1, 0], [0, 1]])
    X_ = np.sqrt(prob) * np.array([[0, 1], [1, 0]])
    return [I_, X_]

def random_unitary(n):
    # draw complex matrix from Ginibre ensemble
    z = np.random.randn(n, n) + 1j * np.random.randn(n, n)
    # QR decompose this complex matrix
    q, r = np.linalg.qr(z)
    # make this decomposition unique
    d = np.diagonal(r)
    l = np.diag(d) / np.abs(d)
    return np.matmul(q, l)

# pick probability
prob = 0.2
# noisy program
p = Program()
p.defgate("DummyGate", random_unitary(2))
p += ("DummyGate", 0)
p.define_noisy_gate("DummyGate", [0], kraus_ops_bit_flip(prob))
ro = p.declare('ro')
p += MEASURE(0, ro)

qc = get_qc('1q-qvm')

num_expts = 1000
num_shots = 1000

p.wrap_in_numshots_loop(num_shots)

results = np.zeros(num_expts)

for i in range(num_expts):
    results_tmp = qc.run(p)
    results[i] = np.count_nonzero(results_tmp) / num_shots

    
plt.figure(figsize=(10, 8))
plt.hist(results, bins=50)
plt.show()

I get the output image centered on the correct value.

Start by reproducing the supplied example and trace the run_and_measure entry point, then compare its execution and measurement path with the explicit run plus MEASURE path. Done means define_noisy_gate produces consistent, unbiased results through run_and_measure, centered on the expected bit-flip probability.