From 5f60c1d9e5ff0fb7df1fb74cecfcca71f02f59f7 Mon Sep 17 00:00:00 2001 From: Asher Yan Date: Sat, 6 Jun 2026 19:37:16 -0400 Subject: [PATCH] Add sampler-backed BackendEstimatorV2 Add EstimatorPub, BackendEstimatorV2, BackendEstimatorJob, and estimator result containers backed by the existing BackendV2 sampler execution path. The implementation extracts Pauli and identity terms from SparseObservable, builds estimator measurement circuits with qubit-wise compatible grouping, derives shots from precision, and post-processes BitArray samples into evs, stds, shots, and target_precision. Include the wrapper fixes needed by the estimator: SparseObservable exposes QkObs array lengths correctly and copies safely, QuantumCircuit supports const-safe copy/measurement inspection and source-circuit compose replay, and hamming_parity/popcount are inline functions for multi-file test builds. --- ...backend-estimator-v2-7575cb4fcb97a8fd.yaml | 11 + src/circuit/quantumcircuit_def.hpp | 18 +- src/primitives/backend_estimator_job.hpp | 288 ++++++++++++++ src/primitives/backend_estimator_v2.hpp | 176 +++++++++ src/primitives/containers/estimator_pub.hpp | 122 ++++++ .../containers/estimator_pub_result.hpp | 86 ++++ .../containers/estimator_result.hpp | 65 ++++ .../detail/estimator_measurement.hpp | 367 ++++++++++++++++++ src/quantum_info/sparse_observable.hpp | 72 +++- src/utils/utils.hpp | 9 +- test/test_estimator.cpp | 190 +++++++++ 11 files changed, 1383 insertions(+), 21 deletions(-) create mode 100644 releasenotes/notes/add-backend-estimator-v2-7575cb4fcb97a8fd.yaml create mode 100644 src/primitives/backend_estimator_job.hpp create mode 100644 src/primitives/backend_estimator_v2.hpp create mode 100644 src/primitives/containers/estimator_pub.hpp create mode 100644 src/primitives/containers/estimator_pub_result.hpp create mode 100644 src/primitives/containers/estimator_result.hpp create mode 100644 src/primitives/detail/estimator_measurement.hpp create mode 100644 test/test_estimator.cpp 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; +}