diff --git a/releasenotes/notes/add-backend-estimator-v2-7575cb4fcb97a8fd.yaml b/releasenotes/notes/add-backend-estimator-v2-7575cb4fcb97a8fd.yaml new file mode 100644 index 0000000..dece9a8 --- /dev/null +++ b/releasenotes/notes/add-backend-estimator-v2-7575cb4fcb97a8fd.yaml @@ -0,0 +1,11 @@ +--- +features: + - | + Add sampler-backed `BackendEstimatorV2`, `EstimatorPub`, + `BackendEstimatorJob`, and estimator result containers. + + The estimator supports real-coefficient Pauli and identity + `SparseObservable` terms for circuits without existing measurements, + derives shots from precision, groups qubit-wise compatible measurement + circuits, and returns result metadata with `target_precision` and `shots`. + Circuits that already contain measurements raise `std::invalid_argument`. diff --git a/src/circuit/quantumcircuit_def.hpp b/src/circuit/quantumcircuit_def.hpp index d8add04..a8673ea 100644 --- a/src/circuit/quantumcircuit_def.hpp +++ b/src/circuit/quantumcircuit_def.hpp @@ -178,10 +178,9 @@ class QuantumCircuit } /// @brief Create a new reference to Quantum Circuit - /// @details Copy constructor of QuantumCircuit does not copy the circuit, - /// but copies shared pointer to Rust's circuit - /// If you want to make a copy of the circuit, - /// please call QuantumCircuit::copy explicitly. + /// @details Copy constructor of QuantumCircuit creates another reference + /// to the same Rust circuit. Call QuantumCircuit::copy explicitly + /// for an independent circuit. /// @param circ a Quantum Circuit to be copied the refrence in the new object QuantumCircuit(const QuantumCircuit &circ) { @@ -260,7 +259,7 @@ class QuantumCircuit /// @brief Copy Quantum Circuit /// @return copied circuit - QuantumCircuit copy(void) + QuantumCircuit copy(void) const { QuantumCircuit copied; copied.rust_circuit_ = std::shared_ptr(qk_circuit_copy(rust_circuit_.get()), qk_circuit_free); @@ -322,6 +321,13 @@ class QuantumCircuit return measure_map_; } + /// @brief get qubits to be measured + /// @return a set of qubits + const std::vector> &get_measure_map(void) const + { + return measure_map_; + } + /// @brief set global phase /// @param phase global phase value void global_phase(const double phase) @@ -1345,7 +1351,7 @@ class QuantumCircuit vclbits[j] = (std::uint32_t)clbits[op->clbits[j]]; } } - QkOperationKind kind = qk_circuit_instruction_kind(rust_circuit_.get(), i); + QkOperationKind kind = qk_circuit_instruction_kind(circ.rust_circuit_.get(), i); if (kind == QkOperationKind_Measure) { qk_circuit_measure(rust_circuit_.get(), vqubits[0], vclbits[0]); } else if (kind == QkOperationKind_Reset) { diff --git a/src/primitives/backend_estimator_job.hpp b/src/primitives/backend_estimator_job.hpp new file mode 100644 index 0000000..104eb64 --- /dev/null +++ b/src/primitives/backend_estimator_job.hpp @@ -0,0 +1,288 @@ +/* +# This code is part of Qiskit. +# +# (C) Copyright IBM 2026. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +*/ + +// job class for BackendEstimator + +#ifndef __qiskitcpp_primitives_backend_estimator_job_hpp__ +#define __qiskitcpp_primitives_backend_estimator_job_hpp__ + +#ifdef _MSC_VER +#define NOMINMAX +#include +#else +#include +#include +#endif + +#include +#include +#include +#include +#include + +#include "primitives/containers/estimator_result.hpp" +#include "primitives/containers/sampler_pub.hpp" +#include "primitives/detail/estimator_measurement.hpp" +#include "providers/job.hpp" + +namespace Qiskit { +namespace primitives { + +/// @class BackendEstimatorJob +/// @brief Job class for Backend Estimator primitive. +class BackendEstimatorJob { + friend class BackendEstimatorV2; +private: + struct PubState { + std::vector evs; + std::vector std_numerators; + uint_t shots; + double target_precision; + }; + + struct JobEntry { + std::shared_ptr job; + uint_t result_index; + estimator_detail::MeasurementTask task; + }; + + std::vector pubs_; + std::vector states_; + std::vector entries_; + bool has_result_ = false; + bool has_submission_error_ = false; + std::string submission_error_; + EstimatorResult cached_result_; + + bool is_final_status(providers::JobStatus status) + { + return status == providers::JobStatus::DONE + || status == providers::JobStatus::CANCELLED + || status == providers::JobStatus::FAILED + || status == providers::JobStatus::ERROR; + } + + std::string status_name(providers::JobStatus status) + { + switch (status) { + case providers::JobStatus::INITIALIZING: + return "INITIALIZING"; + case providers::JobStatus::QUEUED: + return "QUEUED"; + case providers::JobStatus::VALIDATING: + return "VALIDATING"; + case providers::JobStatus::RUNNING: + return "RUNNING"; + case providers::JobStatus::CANCELLED: + return "CANCELLED"; + case providers::JobStatus::DONE: + return "DONE"; + case providers::JobStatus::ERROR: + return "ERROR"; + case providers::JobStatus::FAILED: + return "FAILED"; + } + return "UNKNOWN"; + } + + /// @brief Create a new BackendEstimatorJob. + BackendEstimatorJob(std::vector& pubs, + const std::vector& states, + const std::vector& entries, + const std::string& submission_error = "") + : pubs_(pubs), + states_(states), + entries_(entries), + has_submission_error_(!submission_error.empty()), + submission_error_(submission_error) + { + } + +public: + /// @brief Return the aggregate status of the job. + providers::JobStatus status(void) + { + if (has_submission_error_) + return providers::JobStatus::ERROR; + if (entries_.empty()) + return providers::JobStatus::DONE; + bool has_initializing = false; + bool has_queued = false; + bool has_validating = false; + bool has_running = false; + bool has_cancelled = false; + bool has_failed = false; + bool has_error = false; + for (uint_t i = 0; i < entries_.size(); i++) { + auto st = entries_[i].job->status(); + if (st == providers::JobStatus::ERROR) + has_error = true; + else if (st == providers::JobStatus::FAILED) + has_failed = true; + else if (st == providers::JobStatus::CANCELLED) + has_cancelled = true; + else if (st == providers::JobStatus::RUNNING) + has_running = true; + else if (st == providers::JobStatus::QUEUED) + has_queued = true; + else if (st == providers::JobStatus::VALIDATING) + has_validating = true; + else if (st != providers::JobStatus::DONE) + has_initializing = true; + } + if (has_error) + return providers::JobStatus::ERROR; + if (has_failed) + return providers::JobStatus::FAILED; + if (has_cancelled) + return providers::JobStatus::CANCELLED; + if (has_running) + return providers::JobStatus::RUNNING; + if (has_queued) + return providers::JobStatus::QUEUED; + if (has_validating) + return providers::JobStatus::VALIDATING; + return has_initializing ? providers::JobStatus::INITIALIZING : providers::JobStatus::DONE; + } + + /// @brief Return whether any provider job is actively running. + bool running(void) + { + return status() == providers::JobStatus::RUNNING; + } + + /// @brief Return whether any provider job is queued or validating. + bool queued(void) + { + auto st = status(); + return st == providers::JobStatus::QUEUED || st == providers::JobStatus::VALIDATING; + } + + /// @brief Return whether all provider jobs are done. + bool done(void) + { + return status() == providers::JobStatus::DONE; + } + + /// @brief Return whether any provider job has been cancelled. + bool cancelled(void) + { + return status() == providers::JobStatus::CANCELLED; + } + + /// @brief Return whether every provider job is in a final state. + bool in_final_state(void) + { + return is_final_status(status()); + } + + /// @brief Attempt to cancel the job. + bool cancel(void) + { + return false; + } + + /// @brief Wait for provider jobs and return estimator results. + EstimatorResult result(void) + { + if (has_result_) + return cached_result_; + if (has_submission_error_) { + throw std::runtime_error(submission_error_); + } + + std::vector states = states_; + for (uint_t i = 0; i < entries_.size(); i++) { + providers::JobStatus st; + while (true) { + st = entries_[i].job->status(); + if (is_final_status(st)) + break; +#ifdef _MSC_VER + Sleep(1); +#else + std::this_thread::sleep_for(std::chrono::seconds(1)); +#endif + } + if (st != providers::JobStatus::DONE) { + throw std::runtime_error("BackendEstimatorJob: provider job ended with status " + status_name(st)); + } + if (entries_[i].job->num_results() <= entries_[i].result_index) { + throw std::runtime_error("BackendEstimatorJob: provider result count mismatch"); + } + + circuit::QuantumCircuit result_circuit = entries_[i].task.measured_circuit.copy(); + SamplerPub sampler_pub(result_circuit); + SamplerPubResult sampler_result(sampler_pub); + if (!entries_[i].job->result(entries_[i].result_index, sampler_result)) { + throw std::runtime_error("BackendEstimatorJob: provider result returned false"); + } + + BitArray& bits = sampler_result.data(); + if (bits.num_shots() == 0) { + throw std::runtime_error("BackendEstimatorJob: empty samples for measurement task"); + } + if (bits.num_shots() != entries_[i].task.shots) { + throw std::runtime_error("BackendEstimatorJob: sample count mismatch"); + } + if (bits.num_bits() != entries_[i].task.measured_circuit.num_clbits()) { + throw std::runtime_error("BackendEstimatorJob: sample bit width mismatch"); + } + for (uint_t shot = 0; shot < bits.num_shots(); shot++) { + if (bits[shot].size() != bits.num_bits()) { + throw std::runtime_error("BackendEstimatorJob: sample bit width mismatch"); + } + for (uint_t bit = 0; bit < bits.num_bits(); bit++) { + if (bits[shot][bit] > 1) { + throw std::runtime_error("BackendEstimatorJob: non-binary sample value"); + } + } + } + + uint_t pub_index = entries_[i].task.pub_index; + for (uint_t term_index = 0; term_index < entries_[i].task.group.terms.size(); term_index++) { + const estimator_detail::PauliTerm& term = entries_[i].task.group.terms[term_index]; + auto stats = estimator_detail::pauli_statistics(bits, term); + double variance = stats.variance < 0.0 ? 0.0 : stats.variance; + states[pub_index].evs[term.observable_index] += term.coeff.real() * stats.mean; + states[pub_index].std_numerators[term.observable_index] += std::fabs(term.coeff.real()) * std::sqrt(variance); + } + } + + EstimatorResult result; + for (uint_t pub_index = 0; pub_index < pubs_.size(); pub_index++) { + std::vector stds(states[pub_index].std_numerators.size(), 0.0); + for (uint_t obs_index = 0; obs_index < stds.size(); obs_index++) { + stds[obs_index] = estimator_detail::estimator_std_from_numerator( + states[pub_index].std_numerators[obs_index], + states[pub_index].shots); + } + EstimatorPubResult pub_result(pubs_[pub_index], + states[pub_index].evs, + stds, + states[pub_index].shots, + states[pub_index].target_precision); + result.push_back(pub_result); + } + states_ = states; + cached_result_ = result; + has_result_ = true; + return result; + } +}; + +} // namespace primitives +} // namespace Qiskit + +#endif //__qiskitcpp_primitives_backend_estimator_job_hpp__ diff --git a/src/primitives/backend_estimator_v2.hpp b/src/primitives/backend_estimator_v2.hpp new file mode 100644 index 0000000..d6e94b2 --- /dev/null +++ b/src/primitives/backend_estimator_v2.hpp @@ -0,0 +1,176 @@ +/* +# This code is part of Qiskit. +# +# (C) Copyright IBM 2026. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +*/ + +// estimator implementation for a backend + +#ifndef __qiskitcpp_primitives_backend_estimator_v2_hpp__ +#define __qiskitcpp_primitives_backend_estimator_v2_hpp__ + +#include +#include +#include +#include +#include + +#include "primitives/backend_estimator_job.hpp" +#include "primitives/containers/estimator_pub.hpp" +#include "primitives/containers/sampler_pub.hpp" +#include "primitives/detail/estimator_measurement.hpp" +#include "providers/backend.hpp" + +namespace Qiskit { +namespace primitives { + +/// @class BackendEstimatorV2 +/// @brief Implementation of EstimatorV2 on a sampler-backed backend. +class BackendEstimatorV2 { +protected: + providers::BackendV2& backend_; + double default_precision_; + bool abelian_grouping_; + + std::shared_ptr submit_entries_by_shots( + std::vector& pubs, + const std::vector& states, + std::vector& entries) + { + std::vector submitted(entries.size(), false); + std::vector submitted_entries; + + for (uint_t entry_index = 0; entry_index < entries.size(); entry_index++) { + if (submitted[entry_index]) + continue; + + uint_t shots = entries[entry_index].task.shots; + std::vector sampler_pubs; + std::vector bucket_entries; + + for (uint_t candidate = entry_index; candidate < entries.size(); candidate++) { + if (!submitted[candidate] && entries[candidate].task.shots == shots) { + SamplerPub sampler_pub(entries[candidate].task.measured_circuit); + sampler_pubs.push_back(sampler_pub); + bucket_entries.push_back(candidate); + submitted[candidate] = true; + } + } + + auto provider_job = backend_.run(sampler_pubs, shots); + if (!provider_job) { + return std::shared_ptr( + new BackendEstimatorJob(pubs, states, submitted_entries, "BackendEstimatorV2: backend returned null job")); + } + + for (uint_t result_index = 0; result_index < bucket_entries.size(); result_index++) { + entries[bucket_entries[result_index]].job = provider_job; + entries[bucket_entries[result_index]].result_index = result_index; + submitted_entries.push_back(entries[bucket_entries[result_index]]); + } + } + + return std::shared_ptr(new BackendEstimatorJob(pubs, states, entries)); + } + +public: + /// @brief Create a new BackendEstimatorV2. + /// @param backend backend to run measurement sampler jobs. + /// @param default_precision default target precision. + /// @param abelian_grouping whether to group qubit-wise compatible Pauli terms. + BackendEstimatorV2(providers::BackendV2& backend, + double default_precision = 0.015625, + bool abelian_grouping = true) + : backend_(backend), + default_precision_(default_precision), + abelian_grouping_(abelian_grouping) + { + if (!std::isfinite(default_precision_) || default_precision_ <= 0.0) { + throw std::invalid_argument("BackendEstimatorV2: precision must be greater than 0"); + } + } + + /// @brief return reference to backend object. + const providers::BackendV2& backend(void) const + { + return backend_; + } + + /// @brief Run estimator pubs. + /// @details Input circuits with existing measurements are reported as unsupported. + /// @param pubs An iterable of estimator pub-like objects. + /// @param precision call-level precision. Zero means unset. + /// @return BackendEstimatorJob + /// @throws std::invalid_argument if an input circuit already contains measurements. + std::shared_ptr run(std::vector pubs, + double precision = 0.0) + { + if (!std::isfinite(precision) || precision < 0.0) { + throw std::invalid_argument("BackendEstimatorV2: precision must be non-negative"); + } + + std::vector states(pubs.size()); + std::vector entries; + + for (uint_t pub_index = 0; pub_index < pubs.size(); pub_index++) { + if (!pubs[pub_index].circuit().get_measure_map().empty()) { + throw std::invalid_argument("Estimator measurement: input circuit already contains measurements"); + } + double target_precision = estimator_detail::resolve_precision( + pubs[pub_index].precision(), precision, default_precision_); + uint_t shots = estimator_detail::shots_from_precision(target_precision); + uint_t num_observables = pubs[pub_index].observables().size(); + states[pub_index].evs.resize(num_observables, 0.0); + states[pub_index].std_numerators.resize(num_observables, 0.0); + states[pub_index].shots = shots; + states[pub_index].target_precision = target_precision; + + std::vector terms; + for (uint_t obs_index = 0; obs_index < num_observables; obs_index++) { + if (pubs[pub_index].observables()[obs_index].num_qubits() + != pubs[pub_index].circuit().num_qubits()) { + throw std::invalid_argument("BackendEstimatorV2: observable/circuit width mismatch"); + } + auto extracted = estimator_detail::extract_pauli_terms( + pubs[pub_index].observables()[obs_index], obs_index); + for (uint_t term_index = 0; term_index < extracted.size(); term_index++) { + if (extracted[term_index].factors.empty()) { + states[pub_index].evs[obs_index] += extracted[term_index].coeff.real(); + } else { + terms.push_back(extracted[term_index]); + } + } + } + + auto groups = estimator_detail::group_pauli_terms( + terms, pubs[pub_index].circuit().num_qubits(), abelian_grouping_); + for (uint_t group_index = 0; group_index < groups.size(); group_index++) { + BackendEstimatorJob::JobEntry entry; + entry.result_index = 0; + entry.task.pub_index = pub_index; + entry.task.group_index = group_index; + entry.task.shots = shots; + entry.task.group = groups[group_index]; + entry.task.measured_circuit = estimator_detail::make_measurement_circuit( + pubs[pub_index].circuit(), groups[group_index].basis); + + entries.push_back(entry); + } + } + + return submit_entries_by_shots(pubs, states, entries); + } +}; + +} // namespace primitives +} // namespace Qiskit + +#endif //__qiskitcpp_primitives_backend_estimator_v2_hpp__ diff --git a/src/primitives/containers/estimator_pub.hpp b/src/primitives/containers/estimator_pub.hpp new file mode 100644 index 0000000..a85815c --- /dev/null +++ b/src/primitives/containers/estimator_pub.hpp @@ -0,0 +1,122 @@ +/* +# This code is part of Qiskit. +# +# (C) Copyright IBM 2026. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +*/ + +// estimator pub definition + +#ifndef __qiskitcpp_primitives_estimator_pub_def_hpp__ +#define __qiskitcpp_primitives_estimator_pub_def_hpp__ + +#include +#include +#include + +#include "circuit/quantumcircuit.hpp" +#include "quantum_info/sparse_observable.hpp" + +namespace Qiskit { +namespace primitives { + +/// @class EstimatorPub +/// @brief Estimator Pub(Primitive Unified Bloc). +class EstimatorPub { +protected: + circuit::QuantumCircuit circuit_; + std::vector observables_; + double precision_ = 0.0; + + void validate(void) const + { + if (!std::isfinite(precision_) || precision_ < 0.0) { + throw std::invalid_argument("EstimatorPub: precision must be non-negative"); + } + for (uint_t i = 0; i < observables_.size(); i++) { + if (observables_[i].num_qubits() != circuit_.num_qubits()) { + throw std::invalid_argument("EstimatorPub: observable/circuit width mismatch"); + } + } + } + +public: + /// @brief Create a new EstimatorPub + EstimatorPub() {} + + /// @brief Create a new EstimatorPub + /// @param circuit quantum circuit + /// @param observable observable to estimate + /// @param precision target precision. Zero means unset. + EstimatorPub(const circuit::QuantumCircuit& circuit, + const quantum_info::SparseObservable& observable, + double precision = 0.0) + { + circuit_ = circuit.copy(); + observables_.push_back(observable); + precision_ = precision; + validate(); + } + + /// @brief Create a new EstimatorPub + /// @param circuit quantum circuit + /// @param observables observables to estimate + /// @param precision target precision. Zero means unset. + EstimatorPub(const circuit::QuantumCircuit& circuit, + const std::vector& observables, + double precision = 0.0) + { + circuit_ = circuit.copy(); + observables_ = observables; + precision_ = precision; + validate(); + } + + /// @brief Create a new EstimatorPub as a copy of src. + EstimatorPub(const EstimatorPub& src) + { + circuit_ = src.circuit_.copy(); + observables_ = src.observables_; + precision_ = src.precision_; + } + + EstimatorPub& operator=(const EstimatorPub& src) + { + if (this == &src) + return *this; + circuit_ = src.circuit_.copy(); + observables_ = src.observables_; + precision_ = src.precision_; + return *this; + } + + /// @brief Return a QuantumCircuit for this estimator pub. + const circuit::QuantumCircuit& circuit(void) const + { + return circuit_; + } + + /// @brief Return observables for this estimator pub. + const std::vector& observables(void) const + { + return observables_; + } + + /// @brief Return requested precision. Zero means unset. + double precision(void) const + { + return precision_; + } +}; + +} // namespace primitives +} // namespace Qiskit + +#endif //__qiskitcpp_primitives_estimator_pub_def_hpp__ diff --git a/src/primitives/containers/estimator_pub_result.hpp b/src/primitives/containers/estimator_pub_result.hpp new file mode 100644 index 0000000..826c185 --- /dev/null +++ b/src/primitives/containers/estimator_pub_result.hpp @@ -0,0 +1,86 @@ +/* +# This code is part of Qiskit. +# +# (C) Copyright IBM 2026. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +*/ + +// estimator pub result class + +#ifndef __qiskitcpp_primitives_estimator_pub_result_hpp__ +#define __qiskitcpp_primitives_estimator_pub_result_hpp__ + +#include + +#include "primitives/containers/estimator_pub.hpp" + +namespace Qiskit { +namespace primitives { + +/// @class EstimatorPubResult +/// @brief Result of Estimator Pub(Primitive Unified Bloc). +class EstimatorPubResult { +protected: + EstimatorPub pub_; + std::vector evs_; + std::vector stds_; + uint_t shots_ = 0; + double target_precision_ = 0.0; + +public: + /// @brief Create a new EstimatorPubResult. + EstimatorPubResult() {} + + /// @brief Create a new EstimatorPubResult. + EstimatorPubResult(const EstimatorPub& pub, + const std::vector& evs, + const std::vector& stds, + uint_t shots, + double target_precision) + : pub_(pub), evs_(evs), stds_(stds), shots_(shots), + target_precision_(target_precision) + { + } + + /// @brief Return the pub for this result. + const EstimatorPub& pub(void) const + { + return pub_; + } + + /// @brief Return expectation values. + const std::vector& evs(void) const + { + return evs_; + } + + /// @brief Return standard errors. + const std::vector& stds(void) const + { + return stds_; + } + + /// @brief Return the number of shots used for this pub. + uint_t shots(void) const + { + return shots_; + } + + /// @brief Return resolved target precision for this pub. + double target_precision(void) const + { + return target_precision_; + } +}; + +} // namespace primitives +} // namespace Qiskit + +#endif //__qiskitcpp_primitives_estimator_pub_result_hpp__ diff --git a/src/primitives/containers/estimator_result.hpp b/src/primitives/containers/estimator_result.hpp new file mode 100644 index 0000000..356d323 --- /dev/null +++ b/src/primitives/containers/estimator_result.hpp @@ -0,0 +1,65 @@ +/* +# This code is part of Qiskit. +# +# (C) Copyright IBM 2026. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +*/ + +// estimator result class + +#ifndef __qiskitcpp_primitives_estimator_result_hpp__ +#define __qiskitcpp_primitives_estimator_result_hpp__ + +#include + +#include "primitives/containers/estimator_pub_result.hpp" + +namespace Qiskit { +namespace primitives { + +/// @class EstimatorResult +/// @brief A container for multiple estimator pub results. +class EstimatorResult { +protected: + std::vector pub_results_; + +public: + /// @brief Create a new EstimatorResult. + EstimatorResult() {} + + /// @brief Add an estimator pub result. + void push_back(const EstimatorPubResult& result) + { + pub_results_.push_back(result); + } + + /// @brief Return the number of PUBs in this result. + uint_t size(void) const + { + return pub_results_.size(); + } + + /// @brief Return the pub result. + EstimatorPubResult& operator[](uint_t i) + { + return pub_results_[i]; + } + + /// @brief Return the pub result. + const EstimatorPubResult& operator[](uint_t i) const + { + return pub_results_[i]; + } +}; + +} // namespace primitives +} // namespace Qiskit + +#endif //__qiskitcpp_primitives_estimator_result_hpp__ diff --git a/src/primitives/detail/estimator_measurement.hpp b/src/primitives/detail/estimator_measurement.hpp new file mode 100644 index 0000000..d04d6a5 --- /dev/null +++ b/src/primitives/detail/estimator_measurement.hpp @@ -0,0 +1,367 @@ +/* +# This code is part of Qiskit. +# +# (C) Copyright IBM 2026. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +*/ + +// estimator measurement helpers + +#ifndef __qiskitcpp_primitives_detail_estimator_measurement_hpp__ +#define __qiskitcpp_primitives_detail_estimator_measurement_hpp__ + +#include +#include +#include +#include +#include +#include + +#include "circuit/quantumcircuit.hpp" +#include "primitives/containers/bit_array.hpp" +#include "quantum_info/sparse_observable.hpp" + +namespace Qiskit { +namespace primitives { +namespace estimator_detail { + +constexpr double coefficient_imag_tolerance = 1e-12; + +enum class PauliAxis { + I, + X, + Y, + Z +}; + +struct PauliFactor { + uint_t qubit; + PauliAxis axis; +}; + +struct PauliTerm { + uint_t observable_index; + std::complex coeff; + std::vector factors; +}; + +struct MeasurementBasis { + std::vector axes; +}; + +struct MeasurementGroup { + MeasurementBasis basis; + std::vector terms; +}; + +struct MeasurementTask { + uint_t pub_index; + uint_t group_index; + uint_t shots; + circuit::QuantumCircuit measured_circuit; + MeasurementGroup group; +}; + +struct PauliStats { + double mean; + double variance; +}; + +inline double resolve_precision(double pub_precision, + double run_precision, + double default_precision) +{ + if (pub_precision > 0.0) + return pub_precision; + if (run_precision > 0.0) + return run_precision; + return default_precision; +} + +inline uint_t shots_from_precision(double precision) +{ + if (!std::isfinite(precision) || precision <= 0.0) { + throw std::invalid_argument("BackendEstimatorV2: precision must be greater than 0"); + } + long double p = static_cast(precision); + long double shots = std::ceil(1.0L / (p * p)); + if (!std::isfinite(shots) || shots > static_cast(std::numeric_limits::max())) { + throw std::overflow_error("BackendEstimatorV2: precision maps to shots larger than uint_t"); + } + if (shots < 1.0L) + shots = 1.0L; + return static_cast(shots); +} + +inline PauliAxis axis_from_qk_bit_term(QkBitTerm term) +{ + switch (term) { + case QkBitTerm_X: + return PauliAxis::X; + case QkBitTerm_Y: + return PauliAxis::Y; + case QkBitTerm_Z: + return PauliAxis::Z; + default: + throw std::invalid_argument("Estimator observable: unsupported non-Pauli term"); + } +} + +inline std::vector extract_pauli_terms(const quantum_info::SparseObservable& observable, + uint_t observable_index) +{ + uint_t num_terms = observable.num_terms(); + auto coeffs = observable.coeffs(); + auto bit_terms = observable.bit_terms(); + auto indices = observable.indices(); + auto boundaries = observable.boundaries(); + + if (coeffs.size() != num_terms) { + throw std::runtime_error("Estimator observable: coeffs size mismatch"); + } + if (boundaries.size() != num_terms + 1) { + throw std::runtime_error("Estimator observable: boundaries size mismatch"); + } + if (bit_terms.size() != indices.size()) { + throw std::runtime_error("Estimator observable: bit_terms/indices size mismatch"); + } + if (boundaries.empty() || boundaries[0] != 0 || boundaries.back() != bit_terms.size()) { + throw std::runtime_error("Estimator observable: invalid boundaries"); + } + for (uint_t boundary_index = 1; boundary_index < boundaries.size(); boundary_index++) { + if (boundaries[boundary_index] < boundaries[boundary_index - 1] + || boundaries[boundary_index] > bit_terms.size()) { + throw std::invalid_argument("Estimator observable: invalid boundaries"); + } + } + + std::vector terms; + for (uint_t term_index = 0; term_index < num_terms; term_index++) { + if (std::fabs(coeffs[term_index].imag()) > coefficient_imag_tolerance) { + throw std::invalid_argument("Estimator observable: imaginary coefficient exceeds tolerance"); + } + PauliTerm term; + term.observable_index = observable_index; + term.coeff = coeffs[term_index]; + uint_t begin = boundaries[term_index]; + uint_t end = boundaries[term_index + 1]; + uint_t previous = 0; + for (uint_t pos = begin; pos < end; pos++) { + if (pos > begin && indices[pos] <= previous) { + throw std::runtime_error("Estimator observable: term indices are not strictly increasing"); + } + if (indices[pos] >= observable.num_qubits()) { + throw std::runtime_error("Estimator observable: term index exceeds observable width"); + } + previous = indices[pos]; + PauliFactor factor; + factor.qubit = indices[pos]; + factor.axis = axis_from_qk_bit_term(bit_terms[pos]); + term.factors.push_back(factor); + } + terms.push_back(term); + } + return terms; +} + +inline bool axes_compatible(PauliAxis left, PauliAxis right) +{ + return left == PauliAxis::I || right == PauliAxis::I || left == right; +} + +inline MeasurementBasis term_basis(const PauliTerm& term, uint_t num_qubits) +{ + MeasurementBasis basis; + basis.axes.resize(num_qubits, PauliAxis::I); + for (uint_t i = 0; i < term.factors.size(); i++) { + basis.axes[term.factors[i].qubit] = term.factors[i].axis; + } + return basis; +} + +inline bool term_compatible_with_basis(const PauliTerm& term, const MeasurementBasis& basis) +{ + for (uint_t i = 0; i < term.factors.size(); i++) { + const PauliFactor& factor = term.factors[i]; + if (factor.qubit >= basis.axes.size()) + return false; + if (!axes_compatible(factor.axis, basis.axes[factor.qubit])) + return false; + } + return true; +} + +inline bool terms_compatible(const PauliTerm& left, + const PauliTerm& right, + uint_t num_qubits) +{ + auto basis = term_basis(left, num_qubits); + return term_compatible_with_basis(right, basis); +} + +inline void merge_term_into_basis(const PauliTerm& term, MeasurementBasis& basis) +{ + for (uint_t i = 0; i < term.factors.size(); i++) { + const PauliFactor& factor = term.factors[i]; + if (basis.axes[factor.qubit] == PauliAxis::I) { + basis.axes[factor.qubit] = factor.axis; + } + } +} + +inline bool is_group_valid(const MeasurementGroup& group) +{ + for (uint_t i = 0; i < group.terms.size(); i++) { + if (!term_compatible_with_basis(group.terms[i], group.basis)) + return false; + } + return true; +} + +struct TermOrder { + uint_t index; + uint_t degree; +}; + +inline bool term_order_less(const TermOrder& left, const TermOrder& right) +{ + if (left.degree != right.degree) + return left.degree > right.degree; + return left.index < right.index; +} + +inline std::vector group_pauli_terms(const std::vector& terms, + uint_t num_qubits, + bool abelian_grouping) +{ + std::vector non_identity_terms; + for (uint_t i = 0; i < terms.size(); i++) { + if (!terms[i].factors.empty()) { + non_identity_terms.push_back(terms[i]); + } + } + + if (!abelian_grouping) { + std::vector groups; + for (uint_t i = 0; i < non_identity_terms.size(); i++) { + MeasurementGroup group; + group.basis = term_basis(non_identity_terms[i], num_qubits); + group.terms.push_back(non_identity_terms[i]); + groups.push_back(group); + } + return groups; + } + + std::vector order(non_identity_terms.size()); + for (uint_t i = 0; i < non_identity_terms.size(); i++) { + order[i].index = i; + order[i].degree = 0; + } + for (uint_t i = 0; i < non_identity_terms.size(); i++) { + for (uint_t j = i + 1; j < non_identity_terms.size(); j++) { + if (!terms_compatible(non_identity_terms[i], non_identity_terms[j], num_qubits)) { + order[i].degree++; + order[j].degree++; + } + } + } + std::sort(order.begin(), order.end(), term_order_less); + + std::vector groups; + for (uint_t i = 0; i < order.size(); i++) { + const PauliTerm& term = non_identity_terms[order[i].index]; + bool inserted = false; + for (uint_t group_index = 0; group_index < groups.size(); group_index++) { + if (term_compatible_with_basis(term, groups[group_index].basis)) { + groups[group_index].terms.push_back(term); + merge_term_into_basis(term, groups[group_index].basis); + inserted = true; + break; + } + } + if (!inserted) { + MeasurementGroup group; + group.basis = term_basis(term, num_qubits); + group.terms.push_back(term); + groups.push_back(group); + } + } + return groups; +} + +inline circuit::QuantumCircuit make_measurement_circuit(const circuit::QuantumCircuit& circuit, + const MeasurementBasis& basis) +{ + if (!circuit.get_measure_map().empty()) { + throw std::invalid_argument("Estimator measurement: input circuit already contains measurements"); + } + if (basis.axes.size() != circuit.num_qubits()) { + throw std::invalid_argument("Estimator measurement: basis/circuit width mismatch"); + } + + uint_t num_clbits = std::max(circuit.num_clbits(), circuit.num_qubits()); + circuit::QuantumCircuit measured(circuit.num_qubits(), num_clbits); + circuit::QuantumCircuit source = circuit.copy(); + measured.compose(source); + + for (uint_t qubit = 0; qubit < basis.axes.size(); qubit++) { + if (basis.axes[qubit] == PauliAxis::X) { + measured.h(qubit); + } else if (basis.axes[qubit] == PauliAxis::Y) { + measured.sdg(qubit); + measured.h(qubit); + } + } + for (uint_t qubit = 0; qubit < basis.axes.size(); qubit++) { + if (basis.axes[qubit] != PauliAxis::I) { + measured.measure(qubit, qubit); + } + } + return measured; +} + +inline double pauli_eigenvalue(const BitVector& bits, const PauliTerm& term) +{ + uint_t parity = 0; + for (uint_t i = 0; i < term.factors.size(); i++) { + parity ^= bits[term.factors[i].qubit]; + } + return parity == 0 ? 1.0 : -1.0; +} + +inline PauliStats pauli_statistics(BitArray& bits, const PauliTerm& term) +{ + uint_t shots = bits.num_shots(); + if (shots == 0) { + throw std::runtime_error("BackendEstimatorJob: empty samples for measurement task"); + } + double sum = 0.0; + for (uint_t shot = 0; shot < shots; shot++) { + sum += pauli_eigenvalue(bits[shot], term); + } + PauliStats stats; + stats.mean = sum / static_cast(shots); + stats.variance = 1.0 - stats.mean * stats.mean; + return stats; +} + +inline double estimator_std_from_numerator(double numerator, uint_t shots) +{ + if (shots == 0) { + throw std::runtime_error("BackendEstimatorJob: empty samples for measurement task"); + } + return numerator / std::sqrt(static_cast(shots)); +} + +} // namespace estimator_detail +} // namespace primitives +} // namespace Qiskit + +#endif //__qiskitcpp_primitives_detail_estimator_measurement_hpp__ diff --git a/src/quantum_info/sparse_observable.hpp b/src/quantum_info/sparse_observable.hpp index 0e8895d..203ae87 100644 --- a/src/quantum_info/sparse_observable.hpp +++ b/src/quantum_info/sparse_observable.hpp @@ -19,6 +19,8 @@ #include "qiskit.h" #include "utils/types.hpp" +#include +#include namespace Qiskit { @@ -37,18 +39,55 @@ class SparseObservable { obs_ = nullptr; } - SparseObservable(uint_t num_qubits, std::vector> &coeffs, std::vector &bits, reg_t &indicies, std::vector &boundaries) + SparseObservable(uint_t num_qubits, + std::vector> &coeffs, + std::vector &bits, + reg_t &indices, + std::vector &boundaries) { - std::vector idx32(indicies.size()); - for (int i = 0; i < indicies.size(); i++) { - idx32[i] = (std::uint32_t)indicies[i]; + if (num_qubits > std::numeric_limits::max()) { + throw std::invalid_argument("SparseObservable: number of qubits exceeds QkObs limit"); } - obs_ = qk_obs_new((std::uint32_t)num_qubits, coeffs.size(), bits.size(), (QkComplex64 *)coeffs.data(), bits.data(), idx32.data(), boundaries.data()); + if (bits.size() != indices.size()) { + throw std::invalid_argument("SparseObservable: bits/indices size mismatch"); + } + if (boundaries.size() != coeffs.size() + 1 + || boundaries.empty() + || boundaries[0] != 0 + || boundaries.back() != bits.size()) { + throw std::invalid_argument("SparseObservable: invalid boundaries"); + } + for (uint_t i = 1; i < boundaries.size(); i++) { + if (boundaries[i] < boundaries[i - 1] || boundaries[i] > bits.size()) { + throw std::invalid_argument("SparseObservable: invalid boundaries"); + } + } + for (uint_t term = 0; term < coeffs.size(); term++) { + for (size_t pos = boundaries[term] + 1; pos < boundaries[term + 1]; pos++) { + if (indices[pos] <= indices[pos - 1]) { + throw std::invalid_argument("SparseObservable: term indices are not strictly increasing"); + } + } + } + std::vector idx32(indices.size()); + for (uint_t i = 0; i < indices.size(); i++) { + if (indices[i] >= num_qubits || indices[i] > std::numeric_limits::max()) { + throw std::invalid_argument("SparseObservable: index exceeds observable width"); + } + idx32[i] = (std::uint32_t)indices[i]; + } + obs_ = qk_obs_new((std::uint32_t)num_qubits, + coeffs.size(), + bits.size(), + (QkComplex64 *)coeffs.data(), + bits.data(), + idx32.data(), + boundaries.data()); } SparseObservable(const SparseObservable &other) { - obs_ = qk_obs_copy(other.obs_); + obs_ = other.obs_ ? qk_obs_copy(other.obs_) : nullptr; } ~SparseObservable() @@ -59,6 +98,19 @@ class SparseObservable } } + SparseObservable &operator=(const SparseObservable &other) + { + if (this == &other) { + return *this; + } + if (obs_) { + qk_obs_free(obs_); + obs_ = nullptr; + } + obs_ = other.obs_ ? qk_obs_copy(other.obs_) : nullptr; + return *this; + } + static SparseObservable zero(uint_t num_qubits) { SparseObservable ret; @@ -205,7 +257,7 @@ class SparseObservable } std::vector bit_terms(void) const { - std::vector ret(num_terms()); + std::vector ret(obs_ ? qk_obs_len(obs_) : 0); if (obs_) { auto terms = qk_obs_bit_terms(obs_); @@ -231,7 +283,7 @@ class SparseObservable } reg_t indices(void) const { - reg_t ret(qk_obs_len(obs_)); + reg_t ret(obs_ ? qk_obs_len(obs_) : 0); if (obs_) { auto idx = qk_obs_indices(obs_); @@ -244,7 +296,7 @@ class SparseObservable } reg_t boundaries(void) const { - reg_t ret(qk_obs_len(obs_)); + reg_t ret(obs_ ? qk_obs_num_terms(obs_) + 1 : 0); if (obs_) { auto idx = qk_obs_boundaries(obs_); @@ -310,4 +362,4 @@ class SparseObservable } // namespace quantum_info } // namespace Qiskit -#endif //__qiskitcpp_quantum_info_sparse_observable_hpp__ \ No newline at end of file +#endif //__qiskitcpp_quantum_info_sparse_observable_hpp__ diff --git a/src/utils/utils.hpp b/src/utils/utils.hpp index 6778377..f7f1a61 100644 --- a/src/utils/utils.hpp +++ b/src/utils/utils.hpp @@ -70,15 +70,14 @@ inline bool _naive_parity(uint_t x) { #ifdef INTRINSIC_PARITY -bool (*hamming_parity)(uint_t) = &_intrinsic_parity; -uint_t (*popcount)(uint_t) = &_instrinsic_weight; +inline bool hamming_parity(uint_t x) { return _intrinsic_parity(x); } +inline uint_t popcount(uint_t x) { return _instrinsic_weight(x); } #else -bool (*hamming_parity)(uint_t) = &_naive_parity; -uint_t (*popcount)(uint_t) = &_naive_weight; +inline bool hamming_parity(uint_t x) { return _naive_parity(x); } +inline uint_t popcount(uint_t x) { return _naive_weight(x); } #endif } // namespace Qiskit #endif - diff --git a/test/test_estimator.cpp b/test/test_estimator.cpp new file mode 100644 index 0000000..8c37bf8 --- /dev/null +++ b/test/test_estimator.cpp @@ -0,0 +1,190 @@ +// This code is part of Qiskit. +// +// (C) Copyright IBM 2026. +// +// This code is licensed under the Apache License, Version 2.0. You may +// obtain a copy of this license in the LICENSE.txt file in the root directory +// of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +// +// Any modifications or derivative works of this code must retain this +// copyright notice, and modified files need to carry a notice indicating +// that they have been altered from the originals. + +#include +#include +#include +#include +#include +#include +#include + +#include "common.hpp" + +#include "circuit/quantumcircuit.hpp" +#include "primitives/backend_estimator_v2.hpp" +#include "primitives/containers/estimator_result.hpp" +#include "primitives/detail/estimator_measurement.hpp" +#include "quantum_info/sparse_observable.hpp" +using namespace Qiskit; +using namespace Qiskit::circuit; +using namespace Qiskit::primitives; +using namespace Qiskit::primitives::estimator_detail; +using namespace Qiskit::quantum_info; + +static SparseObservable observable_from_label(const char* text) +{ + std::string label(text); + return SparseObservable::from_label(label); +} +template +static bool raises(Function function) +{ + try { + function(); + } catch (const Exception&) { + return true; + } + return false; +} +static int test_estimator_pub_and_result_containers(void) +{ + QuantumCircuit circuit(1, 0); + auto z = observable_from_label("Z"); + EstimatorPub pub(circuit, z, 0.25); + EstimatorPub assigned; + assigned = pub; + if (assigned.circuit().num_qubits() != 1 || assigned.observables().size() != 1 + || assigned.precision() != 0.25) { + std::cerr << "Unexpected estimator pub contents." << std::endl; + return EqualityError; + } + EstimatorPubResult pub_result(pub, {0.5}, {0.125}, 16, 0.25); + EstimatorResult result; + result.push_back(pub_result); + const EstimatorResult& const_result = result; + if (const_result.size() != 1 || const_result[0].evs()[0] != 0.5 + || const_result[0].stds()[0] != 0.125 || const_result[0].shots() != 16 + || const_result[0].target_precision() != 0.25) { + std::cerr << "Unexpected estimator result contents." << std::endl; + return EqualityError; + } + auto zz = observable_from_label("ZZ"); + if (!raises( + [&] { EstimatorPub(circuit, z, std::numeric_limits::quiet_NaN()); }) + || !raises([&] { EstimatorPub(circuit, zz); })) { + std::cerr << "Estimator pub validation failed." << std::endl; + return EqualityError; + } + return Ok; +} +static int test_precision_and_observable_extraction(void) +{ + if (shots_from_precision(0.015625) != 4096 || shots_from_precision(0.1) != 100 + || std::abs(resolve_precision(0.2, 0.1, 0.015625) - 0.2) > 1e-10 + || std::abs(resolve_precision(0.0, 0.1, 0.015625) - 0.1) > 1e-10) { + std::cerr << "Unexpected precision handling." << std::endl; + return EqualityError; + } + if (!raises([] { shots_from_precision(0.0); })) { + std::cerr << "Zero precision was accepted." << std::endl; + return EqualityError; + } + auto zi = observable_from_label("ZI"); + auto terms = extract_pauli_terms(zi, 0); + if (terms.size() != 1 || terms[0].factors.size() != 1 + || terms[0].factors[0].axis != PauliAxis::Z || terms[0].factors[0].qubit != 1) { + std::cerr << "Unexpected ZI observable extraction." << std::endl; + return EqualityError; + } + std::vector> coeffs{2.0, -0.5}; + std::vector bit_terms{QkBitTerm_Z}; + reg_t indices{0}; + std::vector boundaries{0, 1, 1}; + auto scaled_z_plus_identity = SparseObservable(1, coeffs, bit_terms, indices, boundaries); + auto extracted = extract_pauli_terms(scaled_z_plus_identity, 0); + if (extracted.size() != 2 || extracted[0].factors.size() != 1 + || !extracted[1].factors.empty() || extracted[1].coeff.real() != -0.5) { + std::cerr << "Unexpected multi-term observable extraction." << std::endl; + return EqualityError; + } + auto projector = observable_from_label("+"); + if (!raises([&] { extract_pauli_terms(projector, 0); })) { + std::cerr << "Projector observable was accepted." << std::endl; + return EqualityError; + } + return Ok; +} +static int test_measurement_planning(void) +{ + auto zi = observable_from_label("ZI"); + auto compatible_terms = extract_pauli_terms(zi, 0); + compatible_terms.push_back(PauliTerm{0, 1.0, std::vector{PauliFactor{0, PauliAxis::X}}}); + auto compatible_groups = group_pauli_terms(compatible_terms, 2, true); + if (compatible_groups.size() != 1 || compatible_groups[0].terms.size() != 2 + || !is_group_valid(compatible_groups[0]) + || compatible_groups[0].basis.axes[0] != PauliAxis::X + || compatible_groups[0].basis.axes[1] != PauliAxis::Z) { + std::cerr << "Unexpected compatible Pauli grouping." << std::endl; + return EqualityError; + } + std::vector split_terms{ + PauliTerm{0, 1.0, std::vector{PauliFactor{0, PauliAxis::X}}}, + PauliTerm{0, 1.0, std::vector{PauliFactor{0, PauliAxis::Z}}}}; + auto split_groups = group_pauli_terms(split_terms, 1, true); + if (split_groups.size() != 2 || !is_group_valid(split_groups[0]) + || !is_group_valid(split_groups[1])) { + std::cerr << "Incompatible Pauli terms were grouped together." << std::endl; + return EqualityError; + } + QuantumCircuit circuit(1, 2); + circuit.x(0); + MeasurementBasis y_basis; + y_basis.axes.push_back(PauliAxis::Y); + auto measured = make_measurement_circuit(circuit, y_basis); + if (measured[0].instruction().name() != "x" || measured[1].instruction().name() != "sdg" + || measured[2].instruction().name() != "h" || measured[3].instruction().name() != "measure" + || measured[3].clbits()[0] != 0) { + std::cerr << "Unexpected Y-basis measurement circuit." << std::endl; + return EqualityError; + } + QuantumCircuit premeasured(1, 1); + premeasured.measure(0, 0); + if (!raises([&] { make_measurement_circuit(premeasured, y_basis); })) { + std::cerr << "Pre-measured circuit was accepted." << std::endl; + return EqualityError; + } + return Ok; +} +static int test_pauli_statistics(void) +{ + BitArray samples; + samples.allocate(4, 2); + std::vector bits{reg_t{0, 0}, reg_t{0, 1}, reg_t{0, 1}, reg_t{0, 0}}; + for (uint_t i = 0; i < bits.size(); i++) { + BitVector shot; + shot.from_vector(bits[i]); + samples[i] = shot; + } + PauliTerm zi_term{0, 1.0, std::vector{PauliFactor{1, PauliAxis::Z}}}; + auto stats = pauli_statistics(samples, zi_term); + if (std::abs(stats.mean) > 1e-10 || std::abs(stats.variance - 1.0) > 1e-10 + || std::abs(estimator_std_from_numerator(2.0, 4) - 1.0) > 1e-10) { + std::cerr << "Unexpected Pauli statistics." << std::endl; + return EqualityError; + } + if (!raises([] { estimator_std_from_numerator(1.0, 0); })) { + std::cerr << "Zero-shot standard error was accepted." << std::endl; + return EqualityError; + } + return Ok; +} +int test_estimator(int argc, char** argv) +{ + int num_failed = 0; + num_failed += RUN_TEST(test_estimator_pub_and_result_containers); + num_failed += RUN_TEST(test_precision_and_observable_extraction); + num_failed += RUN_TEST(test_measurement_planning); + num_failed += RUN_TEST(test_pauli_statistics); + std::cerr << "=== Number of failed subtests: " << num_failed << std::endl; + return num_failed; +}