Measure all the observables in a TomographyExperiment.

    :param qc: A QuantumComputer which can run quantum programs
    :param tomo_experiment: A suite of tomographic observables to measure
    :param n_shots: The number of shots to take per ExperimentSetting
    :param progress_callback: If not None, this function is called each time a group of
        settings is run with arguments ``f(i, len(tomo_experiment)`` such that the progress
        is ``i / len(tomo_experiment)``.
    :param active_reset: Whether to actively reset qubits instead of waiting several
        times the coherence length for qubits to decay to |0> naturally. Setting this
        to True is much faster but there is a ~1% error per qubit in the reset operation.
        Thermal noise from "traditional" reset is not routinely characterized but is of the same
        order.
    :param symmetrize_readout: Method used to symmetrize the readout errors, i.e. set
        p(0|1) = p(1|0). For uncorrelated readout errors, this can be achieved by randomly
        selecting between the POVMs {X.D1.X, X.D0.X} and {D0, D1} (where both D0 and D1 are
        diagonal). However, here we currently support exhaustive symmetrization and loop through
        all possible 2^n POVMs {X/I . POVM . X/I}^n, and obtain symmetrization more generally,
        i.e. set p(00|00) = p(01|01) = .. = p(11|11), as well as p(00|01) = p(01|00) etc. If this
        is None, no symmetrization is performed. The exhaustive method can be specified by setting
        this variable to 'exhaustive' (default value). Set to `None` if no symmetrization is
        desired.
    :param calibrate_readout: Method used to calibrate the readout results. Currently, the only
        method supported is normalizing against the operator's expectation value in its +1
        eigenstate, which can be specified by setting this variable to 'plus-eig' (default value).
        The preceding symmetrization and this step together yield a more accurate estimation of the observable. Set to `None` if no calibration is desired.

def measure_observables(qc: QuantumComputer, tomo_experiment: TomographyExperiment,
                        n_shots: int = 10000, progress_callback=None, active_reset=False,
                        symmetrize_readout: str = 'exhaustive',
                        calibrate_readout: str = 'plus-eig'):
    """
    Measure all the observables in a TomographyExperiment.

    :param qc: A QuantumComputer which can run quantum programs
    :param tomo_experiment: A suite of tomographic observables to measure
    :param n_shots: The number of shots to take per ExperimentSetting
    :param progress_callback: If not None, this function is called each time a group of
        settings is run with arguments ``f(i, len(tomo_experiment)`` such that the progress
        is ``i / len(tomo_experiment)``.
    :param active_reset: Whether to actively reset qubits instead of waiting several
        times the coherence length for qubits to decay to |0> naturally. Setting this
        to True is much faster but there is a ~1% error per qubit in the reset operation.
        Thermal noise from "traditional" reset is not routinely characterized but is of the same
        order.
    :param symmetrize_readout: Method used to symmetrize the readout errors, i.e. set
        p(0|1) = p(1|0). For uncorrelated readout errors, this can be achieved by randomly
        selecting between the POVMs {X.D1.X, X.D0.X} and {D0, D1} (where both D0 and D1 are
        diagonal). However, here we currently support exhaustive symmetrization and loop through
        all possible 2^n POVMs {X/I . POVM . X/I}^n, and obtain symmetrization more generally,
        i.e. set p(00|00) = p(01|01) = .. = p(11|11), as well as p(00|01) = p(01|00) etc. If this
        is None, no symmetrization is performed. The exhaustive method can be specified by setting
        this variable to 'exhaustive' (default value). Set to `None` if no symmetrization is
        desired.
    :param calibrate_readout: Method used to calibrate the readout results. Currently, the only
        method supported is normalizing against the operator's expectation value in its +1
        eigenstate, which can be specified by setting this variable to 'plus-eig' (default value).
        The preceding symmetrization and this step together yield a more accurate estimation of the observable. Set to `None` if no calibration is desired.
    """
    # calibration readout only works with symmetrization turned on
    if calibrate_readout is not None and symmetrize_readout is None:
        raise ValueError("Readout calibration only works with readout symmetrization turned on")

    # Outer loop over a collection of grouped settings for which we can simultaneously
    # estimate.
    for i, settings in enumerate(tomo_experiment):

        log.info(f"Collecting bitstrings for the {len(settings)} settings: {settings}")

        # 1.1 Prepare a state according to the amalgam of all setting.in_state
        total_prog = Program()
        if active_reset:
            total_prog += RESET()
        max_weight_in_state = _max_weight_state(setting.in_state for setting in settings)
        for oneq_state in max_weight_in_state.states:
            total_prog += _one_q_state_prep(oneq_state)

        # 1.2 Add in the program
        total_prog += tomo_experiment.program

        # 1.3 Measure the state according to setting.out_operator
        max_weight_out_op = _max_weight_operator(setting.out_operator for setting in settings)
        for qubit, op_str in max_weight_out_op:
            total_prog += _local_pauli_eig_meas(op_str, qubit)

        # 2. Symmetrization
        qubits = max_weight_out_op.get_qubits()

        if symmetrize_readout == 'exhaustive' and len(qubits) > 0:
            bitstrings, d_qub_idx = _exhaustive_symmetrization(qc, qubits, n_shots, total_prog)

        elif symmetrize_readout is None and len(qubits) > 0:
            total_prog_no_symm = total_prog.copy()
            ro = total_prog_no_symm.declare('ro', 'BIT', len(qubits))
            d_qub_idx = {}
            for i, q in enumerate(qubits):
                total_prog_no_symm += MEASURE(q, ro[i])
                # Keep track of qubit-classical register mapping via dict
                d_qub_idx[q] = i
            total_prog_no_symm.wrap_in_numshots_loop(n_shots)
            total_prog_no_symm_native = qc.compiler.quil_to_native_quil(total_prog_no_symm)
            total_prog_no_symm_bin = qc.compiler.native_quil_to_executable(total_prog_no_symm_native)
            bitstrings = qc.run(total_prog_no_symm_bin)

        elif len(qubits) == 0:
            # looks like an identity operation
            pass

        else:
            raise ValueError("Readout symmetrization method must be either 'exhaustive' or None")

        if progress_callback is not None:
            progress_callback(i, len(tomo_experiment))

        # 3. Post-process
        # Inner loop over the grouped settings. They only differ in which qubits' measurements
        # we include in the post-processing. For example, if `settings` is Z1, Z2, Z1Z2 and we
        # measure (n_shots, n_qubits=2) obs_strings then the full operator value involves selecting
        # either the first column, second column, or both and multiplying along the row.
        for setting in settings:
            # 3.1 Get the term's coefficient so we can multiply it in later.
            coeff = complex(setting.out_operator.coefficient)
            if not np.isclose(coeff.imag, 0):
                raise ValueError(f"{setting}'s out_operator has a complex coefficient.")
            coeff = coeff.real

            # 3.2 Special case for measuring the "identity" operator, which doesn't make much
            #     sense but should happen perfectly.
            if is_identity(setting.out_operator):
                yield ExperimentResult(
                    setting=setting,
                    expectation=coeff,
                    std_err=0.0,
                    total_counts=n_shots,
                )
                continue

            # 3.3 Obtain statistics from result of experiment
            obs_mean, obs_var = _stats_from_measurements(bitstrings, d_qub_idx, setting, n_shots, coeff)

            if calibrate_readout == 'plus-eig':
                # 4 Readout calibration
                # 4.1 Obtain calibration program
                calibr_prog = _calibration_program(qc, tomo_experiment, setting)
                # 4.2 Perform symmetrization on the calibration program
                if symmetrize_readout == 'exhaustive':
                    qubs_calibr = setting.out_operator.get_qubits()
                    calibr_shots = n_shots
                    calibr_results, d_calibr_qub_idx = _exhaustive_symmetrization(qc, qubs_calibr, calibr_shots, calibr_prog)

                else:
                    raise ValueError("Readout symmetrization method must be either 'exhaustive' or None")

                # 4.3 Obtain statistics from the measurement process
                obs_calibr_mean, obs_calibr_var = _stats_from_measurements(calibr_results, d_calibr_qub_idx, setting, calibr_shots)
                # 4.3 Calibrate the readout results
                corrected_mean = obs_mean / obs_calibr_mean
                corrected_var = ratio_variance(obs_mean, obs_var, obs_calibr_mean, obs_calibr_var)

                yield ExperimentResult(
                    setting=setting,
                    expectation=corrected_mean.item(),
                    std_err=np.sqrt(corrected_var).item(),
                    total_counts=n_shots,
                    raw_expectation=obs_mean.item(),
                    raw_std_err=np.sqrt(obs_var).item(),
                    calibration_expectation=obs_calibr_mean.item(),
                    calibration_std_err=np.sqrt(obs_calibr_var).item(),
                    calibration_counts=calibr_shots,
                )

            elif calibrate_readout is None:
                # No calibration
                yield ExperimentResult(
                    setting=setting,
                    expectation=obs_mean.item(),
                    std_err=np.sqrt(obs_var).item(),
                    total_counts=n_shots,
                )

            else:
                raise ValueError("Calibration readout method must be either 'plus-eig' or None")