How to use the pyquil.gates.MEASURE function in pyquil

qc = get_qc("9q-generic-qvm")
    compiler = qc.compiler
    qam = qc.qam
    del qc

    q = [4, 5]  # qubits

    # Step 2. Construct your program
    program = Program()
    program += H(q[0])
    program += CNOT(q[0], q[1])

    # Step 2.1. Manage read-out memory
    ro = program.declare("ro", memory_type="BIT", memory_size=2)
    program += MEASURE(q[0], ro[0])
    program += MEASURE(q[1], ro[1])

    # Step 2.2. Run the program in a loop
    program = program.wrap_in_numshots_loop(n_shots)

    # Step 3. Compile and run
    nq_program = compiler.quil_to_native_quil(program)
    executable = compiler.native_quil_to_executable(nq_program)
    bitstrings = qam.load(executable).run().wait().read_memory(region_name="ro")

    # Bincount bitstrings
    basis = np.array([2 ** i for i in range(len(q))])
    ints = np.sum(bitstrings * basis, axis=1)
    print("bincounts", np.bincount(ints))

    # Check parity
    parities = np.sum(bitstrings, axis=1) % 2

for m in qubits:
                program += RX(thetas[var + str(m)], m)

            for m in qubits:
                m_1 = m + 1
                m_sq = m + sq
                skip = (m_1) % sq

                if m_1 < n_qubits:
                    if (m_sq >= n_qubits): program += CNOT(m, m_1)
                    else: program += CNOT(m, m_sq)

                if (m < lim) and (skip != 0): program += CNOT(m, m_1)

            for m in qubits:
                program += MEASURE(m, ro[m])

            print("program instructions from pyquil:\n")
            for instruction in program.instructions:
                print(instruction) 
                
            return program

print("\nMinimum energy =", round(result["fun"], 4))
print("Best angle =", round(result["x"][0], 4))

# =======================================
# VQE on a noisy simulator: Pauli Channel
# =======================================
print()
print(" VQE on a noisy simulator: Pauli channel ".center(80, "="))
print()

# Create a noise model which has a 10% chance of each gate at each timestep
pauli_channel = [0.1, 0.1, 0.1]
noisy_qvm = api.QVMConnection(gate_noise=pauli_channel)

# Check that the simulator is indeed noisy
p = Program(X(0), MEASURE(0, 0))
res = noisy_qvm.run(p, [0], 10)
print("Outcome of NOT and MEASURE circuit on noisy simulator:")
print(res)

# Update the minimizer in VQE to start with a larger initial simplex
vqe_inst.minimizer_kwargs = {"method": "Nelder-mead", 
                             "options": 
                                 {"initial_simplex": np.array([[0.0], [0.05]]), 
                                  "xatol": 1.0e-2}
                                 }

# Loop over a range of angles and plot expectation with sampling 
sampled = [vqe_inst.expectation(small_ansatz([angle]), hamiltonian, 1000, noisy_qvm)
           for angle in angle_range]

# Plot the sampled expectation

and accumulating the resulting bitarrays.

    :param qc: a quantum computer object on which to run each program
    :param programs: a list of programs to run
    :param meas_qubits: groups of qubits to measure for each program
    :param num_shots: the number of shots to run for each program
    :return: a len(programs) long list of num_shots by num_meas_qubits bit arrays of results for
        each program.
    """
    results = []
    for program in programs:
        # copy the program so the original is not mutated
        prog = program.copy()
        ro = prog.declare("ro", "BIT", len(meas_qubits))
        for idx, q in enumerate(meas_qubits):
            prog += MEASURE(q, ro[idx])

        prog.wrap_in_numshots_loop(num_shots)
        prog = qc.compiler.quil_to_native_quil(prog)
        exe = qc.compiler.native_quil_to_executable(prog)
        shots = qc.run(exe)
        results.append(shots)
    return results

Runs Grover's Algorithm to find the bitstring that is designated by ``bistring_map``.

        In particular, this will prepare an initial state in the uniform superposition over all bit-
        strings, an then use Grover's Algorithm to pick out the desired bitstring.

        :param qc: the connection to the Rigetti cloud to run pyQuil programs.
        :param bitstring_map: a mapping from bitstrings to the phases that the oracle should impart
            on them. If the oracle should "look" for a bitstring, it should have a ``-1``, otherwise
            it should have a ``1``.
        :return: Returns the bitstring resulting from measurement after Grover's Algorithm.
        """

        self._init_attr(bitstring_map)

        ro = self.grover_circuit.declare('ro', 'BIT', len(self.qubits))
        self.grover_circuit += [MEASURE(qubit, ro[idx]) for idx, qubit in enumerate(self.qubits)]
        executable = qc.compile(self.grover_circuit)
        sampled_bitstring = qc.run(executable)

        return "".join([str(bit) for bit in sampled_bitstring[0]])

In contrast to :py:class:`QVMConnection.run_and_measure`, this method simulates
            noise correctly for noisy QVMs. However, this method is slower for ``trials > 1``.
            For faster noise-free simulation, consider
            :py:class:`WavefunctionSimulator.run_and_measure`.

        :param program: The state preparation program to run and then measure.
        :param trials: The number of times to run the program.
        :return: A dictionary keyed by qubit index where the corresponding value is a 1D array of
            measured bits.
        """
        program = program.copy()
        validate_protoquil(program)
        ro = program.declare('ro', 'BIT', len(self.qubits()))
        for i, q in enumerate(self.qubits()):
            program.inst(MEASURE(q, ro[i]))
        program.wrap_in_numshots_loop(trials)
        executable = self.compile(program)
        bitstring_array = self.run(executable=executable)
        bitstring_dict = {}
        for i, q in enumerate(self.qubits()):
            bitstring_dict[q] = bitstring_array[:, i]
        return bitstring_dict

Given a wavefunctions, this calculates the expectation value of the Zi
    operator where i ranges over all the qubits given in marked_qubits.

    :param pyquil_program: pyQuil program generating some state
    :param marked_qubits: The qubits within the support of the Z pauli
                          operator whose expectation value is being calculated
    :param qc: A QuantumComputer object.
    :param samples: Number of bitstrings collected to calculate expectation
                    from sampling.
    :returns: The expectation value as a float.
    """
    program = Program()
    ro = program.declare('ro', 'BIT', max(marked_qubits) + 1)
    program += pyquil_program
    program += [MEASURE(qubit, r) for qubit, r in zip(list(range(max(marked_qubits) + 1)), ro)]
    program.wrap_in_numshots_loop(samples)
    executable = qc.compile(program)
    bitstring_samples = qc.run(executable)
    bitstring_tuples = list(map(tuple, bitstring_samples))

    freq = Counter(bitstring_tuples)

    # perform weighted average
    expectation = 0
    for bitstring, count in freq.items():
        bitstring_int = int("".join([str(x) for x in bitstring[::-1]]), 2)
        if parity_even_p(bitstring_int, marked_qubits):
            expectation += float(count) / samples
        else:
            expectation -= float(count) / samples
    return expectation

def _wrap_program(self, program):
        # the actions select gates. but a pyquil program needs a bit more
        # namely, declaration of classical memory for readout, and suitable
        # measurement instructions
        ro = program.declare('ro', 'BIT', self.num_qubits)
        for q in range(self.num_qubits):
            program.inst(MEASURE(q, ro[q]))
        program.wrap_in_numshots_loop(NUM_SHOTS)
        return program

pauli_terms = pauli_terms.terms

    # check if each term commutes with everything
    if commutation_check:
        if len(commuting_sets(sum(pauli_terms))) != 1:
            raise CommutationError("Not all terms commute in the expected way")

    program = copy.deepcopy(program)
    pauli_for_rotations = PauliTerm.from_list(
        [(value, key) for key, value in basis_transform_dict.items()])

    program += get_rotation_program(pauli_for_rotations)

    qubits = sorted(list(basis_transform_dict.keys()))
    for qubit in qubits:
        program.inst(MEASURE(qubit, qubit))

    coeff_vec = np.array(
        list(map(lambda x: x.coefficient, pauli_terms))).reshape((-1, 1))

    # upper bound on samples given by IV of arXiv:1801.03524
    num_sample_ubound = int(np.ceil(np.sum(np.abs(coeff_vec))**2 / variance_bound))
    results = None
    sample_variance = np.infty
    number_of_samples = 0
    while (sample_variance > variance_bound and
           number_of_samples < num_sample_ubound):
        ro = program.declare('ro', 'BIT', len(qubits))
        program += [MEASURE(qubit, ro_reg) for qubit, ro_reg in zip(qubits, ro)]
        program.wrap_in_numshots_loop(min(10000, num_sample_ubound))
        executable = quantum_resource.compiler.native_quil_to_executable(program)
        tresults = quantum_resource.run(executable)