diff --git a/README.md b/README.md index a12fbe0..801d80e 100644 --- a/README.md +++ b/README.md @@ -92,6 +92,9 @@ pyqpanda-algorithm 是由本源量子（Origin Quantum）开发的量子算法 - **SparseAmp（稀疏幅度编码）** 高效将稀疏向量编码为量子态，减少量子资源消耗，适用于数据预处理。 +- **QShadow（可观测量感知的经典阴影）** + 通过局域偏置 Pauli 测量从同一批样本估计多个可观测量，并提供同时置信区间与测量预算控制。 + ------ ## 安装 diff --git a/README_EN.md b/README_EN.md index d74c953..4c0ed36 100644 --- a/README_EN.md +++ b/README_EN.md @@ -59,6 +59,8 @@ Provides underlying tools for quantum computing workflows. Implements numerical comparison or threshold judgment to build quantum decision logic. - **SparseAmp (Sparse Amplitude Encoding)** Efficiently encodes sparse vectors into quantum states, reducing quantum resource consumption (for data preprocessing). +- **QShadow (Observable-Aware Classical Shadows)** + Estimates many observables from locally biased Pauli measurements with simultaneous confidence intervals and explicit shot-budget control. ------ diff --git a/Tutorials/source/QShadow.rst b/Tutorials/source/QShadow.rst new file mode 100644 index 0000000..f2a8a8c --- /dev/null +++ b/Tutorials/source/QShadow.rst @@ -0,0 +1,96 @@ +QShadow：可观测量感知的经典阴影 +==================================== + +QShadow 从同一批局域 Pauli 测量中估计多个可观测量。它适合变分算法能量读出、 +相关函数扫描和小规模真机实验的测量后处理。经典阴影与局域偏置方法是已有研究； +本模块的 PyQPanda3 集成、带下限规划器、置信预算流程和测试为独立实现。 + +基本工作流 +---------- + +.. code-block:: python + + from pyqpanda_alg.QShadow import ( + PauliObservable, + optimize_basis_probabilities, + sample_measurement_plan, + simulate_pauli_measurements, + estimate_observables, + ) + + targets = [ + PauliObservable.from_terms( + {"ZZ": -1.0, "XI": 0.2, "IX": 0.2}, name="energy" + ), + PauliObservable.from_terms({"ZI": 1.0}, name="magnetisation"), + ] + result = optimize_basis_probabilities(targets, probability_floor=0.03) + plan = sample_measurement_plan(result.probabilities, 5000, seed=7) + data = simulate_pauli_measurements([1, 0, 0, 0], plan, seed=8) + estimates = estimate_observables(data, targets, family_confidence=0.95) + +约定与无偏性 +------------ + +Pauli 字符串和数据集统一使用 q0-first：``"XI"`` 表示 X 作用于 q0。对目标 Pauli +项，只有测量基覆盖其非恒等支撑时，该 shot 才产生非零贡献；贡献除以对应测量基 +概率的乘积。数据集保存每个 shot 的完整 X/Y/Z 概率，因此不同规划轮次可以无偏 +合并。 + +规划器优化的是系数加权的对角二阶矩代理，并通过 ``probability_floor`` 防止某个 +基完全不可观测及逆倾向权重发散。代理损失下降不是任意未知态上实际方差下降的 +保证。 + +离线重复基准 +-------------- + +``example/QShadow/qshadow_planning_benchmark.py`` 使用固定状态、可观测量和种子， +对均匀规划与 QShadow 规划进行 1000 轮、每轮 2000 shots 的对照。规划损失从 +``24.54`` 降至 ``5.964549``，聚合 RMSE 降低 ``60.20%``，两个可观测量的平均 +区间半宽分别降低 ``62.35%`` 和 ``53.42%``。在另一个 100 轮顺序预算实验中， +达到 ``0.22`` 目标半宽所需平均 shots 从 ``7284.9`` 降至 ``1892.8``，减少 +``74.02%``，且没有硬上限违规。 + +两种方案在该配置中的单项和 family-wise 经验覆盖率均为 ``1.000``，说明区间较 +保守，不表示精确 95% 校准。结果只适用于该固定基准，不构成所有状态与可观测量上 +普遍优势的声明。 + +置信预算 +-------- + +``ShadowBudgetController`` 用 Bonferroni 同时经验 Bernstein 半宽决定下一批 shots， +同时遵守 ``max_shots``。它只返回建议，不自动提交任务。远程执行还可通过 +``run_measurement_plan(..., max_jobs=N)`` 在任何任务提交前限制唯一测量线路数量。 + +真机结果必须同时报告设备/后端、时间、shots、任务 ID、原始 counts、采样概率、 +统计区间及理想模拟器偏差。当前实现不含读出误差缓解，因此经验 Bernstein 区间 +不覆盖硬件系统偏差。 + +WK_C180 实测证据 +---------------- + +2026-07-26 的最小验证在 ``WK_C180`` 上执行 X/Y/Z 各 512 shots，总计恰好 +3 个任务、1536 shots，且没有自动重试。真实 counts 经 q0-first 转换后进入同一个 +QShadow 数据集，按逐 shot 的 ``(1/3, 1/3, 1/3)`` 概率同时估计 X/Y/Z，并生成 +family-wise 经验 Bernstein 区间。机器可读证据位于 +``example/QShadow/qshadow_wukong_validation_20260726.json``，配套解释位于同目录 +``qshadow_wukong_validation_20260726.md``。 + +这次运行验证了从真实后端 counts 到同时估计、置信区间和 3-job/1536-shot 预算的 +完整工作流。硬件 X/Y/Z 区间均未覆盖理想值，Bloch 向量范数为 +``1.539826 > 1``；本次平台配置包含 ``is_amend=True``，QShadow 未执行读出误差 +缓解，因此不把结果解释为硬件精度成功。由于 X/Y/Z 目标完全对称并采用均匀计划， +本实验也不比较偏置规划器与均匀规划的真机优势；规划器行为由 CPUQVM 示例和测试 +覆盖。 + +引用 +---- + +* Huang, Kueng, Preskill, *Predicting many properties of a quantum system from + very few measurements*, Nature Physics 16 (2020), arXiv:2002.08953. +* Hadfield, Bravyi, Raymond, Mezzacapo, *Measurements of Quantum Hamiltonians + with Locally-Biased Classical Shadows*, arXiv:2006.15788. +* Hadfield, *Adaptive Pauli Shadows for Energy Estimation*, arXiv:2105.12207. + +完整 API、云后端安全示例、限制和原创性边界见 +``pyqpanda_alg/QShadow/README.md``。 diff --git a/Tutorials/source/index.rst b/Tutorials/source/index.rst index 0866854..32a3582 100644 --- a/Tutorials/source/index.rst +++ b/Tutorials/source/index.rst @@ -24,6 +24,12 @@ Overall, it provides a standardized set of tools for developers, allowing them t Changelog +.. toctree:: + :maxdepth: 2 + :caption: QShadow Tutorial + + QShadow + .. toctree:: :caption: API Reference :maxdepth: 2 @@ -42,3 +48,4 @@ Overall, it provides a standardized set of tools for developers, allowing them t autoapi/pyqpanda_alg/QSVD/index autoapi/pyqpanda_alg/QSVR/index autoapi/pyqpanda_alg/QUBO/index + autoapi/pyqpanda_alg/QShadow/index diff --git a/pyqpanda-algorithm/example/QShadow/qshadow_observable_estimation.py b/pyqpanda-algorithm/example/QShadow/qshadow_observable_estimation.py new file mode 100644 index 0000000..0e61e97 --- /dev/null +++ b/pyqpanda-algorithm/example/QShadow/qshadow_observable_estimation.py @@ -0,0 +1,106 @@ +"""QShadow end-to-end example: planning, CPUQVM sampling, and confidence budget.""" + +from __future__ import annotations + +from pyqpanda3.core import CNOT, CPUQVM, QCircuit, QProg, RY + +from pyqpanda_alg.QShadow import ( + PauliObservable, + PyQPandaRunner, + ShadowBudgetController, + estimate_observables, + optimize_basis_probabilities, + planning_loss, + run_measurement_plan, + sample_measurement_plan, + uniform_probabilities, +) + + +def prepare(qubits): + """A reproducible three-qubit entangled state-preparation circuit.""" + + return ( + QCircuit() + << RY(qubits[0], 0.7) + << RY(qubits[1], -0.4) + << RY(qubits[2], 0.3) + << CNOT(qubits[0], qubits[1]) + << CNOT(qubits[1], qubits[2]) + ) + + +def exact_statevector(): + program = QProg(3) + program << prepare(program.qubits()) + machine = CPUQVM() + machine.run(program, shots=1) + return machine.result().get_state_vector() + + +def main(): + observables = [ + PauliObservable.from_terms( + { + "ZZI": -1.0, + "IZZ": -0.8, + "XII": 0.20, + "IXI": 0.15, + "IIX": 0.10, + }, + name="structured_energy", + ), + PauliObservable.from_terms( + {"ZIZ": 1.0, "XXI": 0.25}, name="long_range_order" + ), + ] + + uniform = uniform_probabilities(3) + optimisation = optimize_basis_probabilities( + observables, probability_floor=0.04 + ) + print("uniform planning loss :", planning_loss(uniform, observables)) + print("optimised planning loss:", optimisation.final_loss) + print("basis probabilities [X,Y,Z]:\n", optimisation.probabilities) + + plan = sample_measurement_plan( + optimisation.probabilities, shots=4_000, seed=20260726 + ) + print("shots / unique backend jobs:", plan.shots, len(plan.grouped_shots)) + + runner = PyQPandaRunner(prepare, n_qubits=3) + dataset = run_measurement_plan(plan, runner, max_jobs=27) + estimates = estimate_observables( + dataset, observables, family_confidence=0.95 + ) + state = exact_statevector() + for observable in observables: + estimate = estimates[observable.name] + exact = observable.exact_expectation(state) + print( + f"{observable.name}: estimate={estimate.value:.5f}, " + f"exact={exact:.5f}, CI=[{estimate.lower:.5f}, {estimate.upper:.5f}], " + f"min_matches={estimate.minimum_term_matches}" + ) + + controller = ShadowBudgetController( + observables, + target_half_width=0.25, + family_confidence=0.95, + min_shots=1_000, + max_shots=8_000, + batch_shots=500, + ) + report = controller.evaluate(dataset) + print( + "budget:", + report.reason, + "next_batch_shots=", + report.next_batch_shots, + "hard_max=", + report.max_shots, + ) + + +if __name__ == "__main__": + main() diff --git a/pyqpanda-algorithm/example/QShadow/qshadow_planning_benchmark.py b/pyqpanda-algorithm/example/QShadow/qshadow_planning_benchmark.py new file mode 100644 index 0000000..fbe7bef --- /dev/null +++ b/pyqpanda-algorithm/example/QShadow/qshadow_planning_benchmark.py @@ -0,0 +1,378 @@ +"""Reproducible offline benchmark for QShadow planning and shot budgets. + +This script compares uniform and observable-aware local-Pauli plans on the same +state, observables, shot counts, confidence level, and fixed seed blocks. It +uses CPUQVM once to obtain the exact statevector and NumPy sampling afterwards; +it never connects to a cloud backend. +""" + +from __future__ import annotations + +import argparse +import json +from collections.abc import Sequence +from pathlib import Path +from typing import Any + +import numpy as np +from pyqpanda3.core import CNOT, CPUQVM, QCircuit, QProg, RY + +from pyqpanda_alg.QShadow import ( + PauliObservable, + ShadowBudgetController, + ShadowDataset, + estimate_observables, + optimize_basis_probabilities, + planning_loss, + sample_measurement_plan, + simulate_pauli_measurements, + uniform_probabilities, +) + + +def prepare(qubits): + """Prepare the same three-qubit state used by the QShadow example.""" + + return ( + QCircuit() + << RY(qubits[0], 0.7) + << RY(qubits[1], -0.4) + << RY(qubits[2], 0.3) + << CNOT(qubits[0], qubits[1]) + << CNOT(qubits[1], qubits[2]) + ) + + +def exact_statevector() -> list[complex]: + """Return the CPUQVM statevector before measurement.""" + + program = QProg(3) + program << prepare(program.qubits()) + machine = CPUQVM() + machine.run(program, shots=1) + return [complex(value) for value in machine.result().get_state_vector()] + + +def benchmark_observables() -> list[PauliObservable]: + """Return the fixed observable family used by this benchmark.""" + + return [ + PauliObservable.from_terms( + { + "ZZI": -1.0, + "IZZ": -0.8, + "XII": 0.20, + "IXI": 0.15, + "IIX": 0.10, + }, + name="structured_energy", + ), + PauliObservable.from_terms( + {"ZIZ": 1.0, "XXI": 0.25}, name="long_range_order" + ), + ] + + +def _fixed_shot_statistics( + *, + state: Sequence[complex], + observables: list[PauliObservable], + probabilities: dict[str, list[list[float]]], + trials: int, + shots: int, + confidence: float, + seed: int, +) -> dict[str, Any]: + exact = np.asarray( + [observable.exact_expectation(state) for observable in observables] + ) + results: dict[str, Any] = {} + for method_index, (method, method_probabilities) in enumerate( + probabilities.items() + ): + values: list[list[float]] = [] + half_widths: list[list[float]] = [] + covered: list[list[bool]] = [] + for trial in range(trials): + offset = method_index * 200_000 + trial + plan = sample_measurement_plan( + method_probabilities, shots, seed=seed + offset + ) + dataset = simulate_pauli_measurements( + state, plan, seed=seed + 500_000 + offset + ) + estimates = estimate_observables( + dataset, observables, family_confidence=confidence + ) + values.append( + [estimates[observable.name].value for observable in observables] + ) + half_widths.append( + [ + estimates[observable.name].half_width + for observable in observables + ] + ) + covered.append( + [ + estimates[observable.name].lower + <= exact_value + <= estimates[observable.name].upper + for observable, exact_value in zip(observables, exact) + ] + ) + + value_array = np.asarray(values) + width_array = np.asarray(half_widths) + covered_array = np.asarray(covered, dtype=bool) + errors = value_array - exact + per_observable: dict[str, Any] = {} + for index, observable in enumerate(observables): + bias = float(np.mean(errors[:, index])) + bias_standard_error = float( + np.std(errors[:, index], ddof=1) / np.sqrt(trials) + ) + per_observable[observable.name] = { + "bias": bias, + "bias_standard_error": bias_standard_error, + "bias_z_score": bias / bias_standard_error, + "rmse": float(np.sqrt(np.mean(errors[:, index] ** 2))), + "mean_absolute_error": float(np.mean(np.abs(errors[:, index]))), + "mean_half_width": float(np.mean(width_array[:, index])), + "empirical_coverage": float(np.mean(covered_array[:, index])), + } + results[method] = { + "per_observable": per_observable, + "aggregate_rmse": float(np.sqrt(np.mean(errors**2))), + "family_wise_empirical_coverage": float( + np.mean(np.all(covered_array, axis=1)) + ), + } + return {"exact_values": dict(zip((o.name for o in observables), exact)), **results} + + +def _sequential_budget_statistics( + *, + state: Sequence[complex], + observables: list[PauliObservable], + probabilities: dict[str, list[list[float]]], + trials: int, + target_half_width: float, + confidence: float, + min_shots: int, + max_shots: int, + batch_shots: int, + seed: int, +) -> dict[str, Any]: + exact = { + observable.name: observable.exact_expectation(state) + for observable in observables + } + results: dict[str, Any] = {} + for method_index, (method, method_probabilities) in enumerate( + probabilities.items() + ): + final_shots: list[int] = [] + rounds: list[int] = [] + completed: list[bool] = [] + family_covered: list[bool] = [] + hard_max_violations = 0 + width_rule_violations = 0 + for trial in range(trials): + controller = ShadowBudgetController( + observables, + target_half_width=target_half_width, + family_confidence=confidence, + min_shots=min_shots, + max_shots=max_shots, + batch_shots=batch_shots, + ) + datasets: list[ShadowDataset] = [] + next_shots = min_shots + step = 0 + while True: + offset = method_index * 500_000 + trial * 50 + step + plan = sample_measurement_plan( + method_probabilities, + next_shots, + seed=seed + 2_000_000 + offset, + ) + datasets.append( + simulate_pauli_measurements( + state, + plan, + seed=seed + 3_000_000 + offset, + ) + ) + dataset = ShadowDataset.concatenate(*datasets) + report = controller.evaluate(dataset) + step += 1 + if dataset.shots > max_shots: + hard_max_violations += 1 + if report.complete: + if any( + estimate.half_width > target_half_width + 1e-12 + for estimate in report.estimates.values() + ): + width_rule_violations += 1 + break + if report.exhausted: + break + next_shots = report.next_batch_shots + if step > 20: + raise RuntimeError("budget controller did not terminate") + + final_shots.append(dataset.shots) + rounds.append(step) + completed.append(report.complete) + family_covered.append( + all( + estimate.lower <= exact[name] <= estimate.upper + for name, estimate in report.estimates.items() + ) + ) + + shot_array = np.asarray(final_shots) + results[method] = { + "target_reached_rate": float(np.mean(completed)), + "family_wise_empirical_coverage": float(np.mean(family_covered)), + "hard_max_violations": hard_max_violations, + "width_rule_violations": width_rule_violations, + "final_shots": { + "mean": float(np.mean(shot_array)), + "median": float(np.median(shot_array)), + "p10": float(np.quantile(shot_array, 0.1)), + "p90": float(np.quantile(shot_array, 0.9)), + "minimum": int(np.min(shot_array)), + "maximum": int(np.max(shot_array)), + }, + "mean_rounds": float(np.mean(rounds)), + } + return results + + +def run_benchmark( + *, + trials: int = 200, + shots: int = 2_000, + budget_trials: int = 100, + target_half_width: float = 0.22, + confidence: float = 0.95, + min_shots: int = 500, + max_shots: int = 12_000, + batch_shots: int = 250, + seed: int = 2_000_000, +) -> dict[str, Any]: + """Run the fixed-shot and sequential-budget benchmark.""" + + if trials <= 1 or budget_trials <= 0 or shots <= 0: + raise ValueError("trials must exceed one; budget_trials and shots must be positive") + state = exact_statevector() + observables = benchmark_observables() + uniform = uniform_probabilities(3).tolist() + optimized = optimize_basis_probabilities( + observables, probability_floor=0.04 + ).probabilities.tolist() + probabilities = {"uniform": uniform, "optimized": optimized} + fixed = _fixed_shot_statistics( + state=state, + observables=observables, + probabilities=probabilities, + trials=trials, + shots=shots, + confidence=confidence, + seed=seed, + ) + budget = _sequential_budget_statistics( + state=state, + observables=observables, + probabilities=probabilities, + trials=budget_trials, + target_half_width=target_half_width, + confidence=confidence, + min_shots=min_shots, + max_shots=max_shots, + batch_shots=batch_shots, + seed=seed, + ) + + uniform_fixed = fixed["uniform"] + optimized_fixed = fixed["optimized"] + rmse_reduction = { + observable.name: 1.0 + - optimized_fixed["per_observable"][observable.name]["rmse"] + / uniform_fixed["per_observable"][observable.name]["rmse"] + for observable in observables + } + width_reduction = { + observable.name: 1.0 + - optimized_fixed["per_observable"][observable.name]["mean_half_width"] + / uniform_fixed["per_observable"][observable.name]["mean_half_width"] + for observable in observables + } + shot_reduction = 1.0 - ( + budget["optimized"]["final_shots"]["mean"] + / budget["uniform"]["final_shots"]["mean"] + ) + return { + "configuration": { + "fixed_shot_trials": trials, + "shots_per_fixed_trial": shots, + "budget_trials": budget_trials, + "target_half_width": target_half_width, + "family_confidence": confidence, + "min_shots": min_shots, + "max_shots": max_shots, + "batch_shots": batch_shots, + "seed": seed, + }, + "planning_loss": { + "uniform": planning_loss(uniform, observables), + "optimized": planning_loss(optimized, observables), + }, + "fixed_shot_statistics": fixed, + "sequential_budget_statistics": budget, + "optimized_reduction_fraction": { + "aggregate_rmse": 1.0 + - optimized_fixed["aggregate_rmse"] + / uniform_fixed["aggregate_rmse"], + "rmse_by_observable": rmse_reduction, + "mean_half_width_by_observable": width_reduction, + "mean_final_shots_to_target": shot_reduction, + }, + "interpretation": ( + "Observed coverage is an empirical diagnostic, not proof of exact " + "95% calibration. The comparison is specific to this fixed state, " + "observable family, probability floor, and seed schedule." + ), + } + + +def _parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--trials", type=int, default=200) + parser.add_argument("--shots", type=int, default=2_000) + parser.add_argument("--budget-trials", type=int, default=100) + parser.add_argument("--target-half-width", type=float, default=0.22) + parser.add_argument("--seed", type=int, default=2_000_000) + parser.add_argument("--output", type=Path) + return parser.parse_args() + + +def main() -> None: + args = _parse_args() + results = run_benchmark( + trials=args.trials, + shots=args.shots, + budget_trials=args.budget_trials, + target_half_width=args.target_half_width, + seed=args.seed, + ) + payload = json.dumps(results, indent=2, ensure_ascii=False) + if args.output is not None: + args.output.write_text(payload + "\n", encoding="utf-8") + print(payload) + + +if __name__ == "__main__": + main() diff --git a/pyqpanda-algorithm/example/QShadow/qshadow_wukong_validation_20260726.json b/pyqpanda-algorithm/example/QShadow/qshadow_wukong_validation_20260726.json new file mode 100644 index 0000000..726fcae --- /dev/null +++ b/pyqpanda-algorithm/example/QShadow/qshadow_wukong_validation_20260726.json @@ -0,0 +1,272 @@ +{ + "schema_version": 1, + "run_id": "qshadow-wukong-xyz-20260726", + "created_at_beijing": "2026-07-26T12:31:43+08:00", + "backend": "WK_C180", + "service_url": "https://qpanda3-runtime.qpanda.cn", + "sdk_versions": { + "qpanda3_runtime": "1.0.0", + "pyqpanda3": "0.3.4" + }, + "budget": { + "max_jobs": 3, + "shots_per_job": 512, + "max_total_shots": 1536, + "automatic_retry": false + }, + "state_preparation": { + "initial_state": "|0>", + "gates": [ + "RY(0.9) q[0]", + "RZ(0.7) q[0]" + ] + }, + "measurement_plan": { + "q0_first": true, + "basis_order": [ + "X", + "Y", + "Z" + ], + "basis_probabilities_xyz": [ + 0.3333333333333333, + 0.3333333333333333, + 0.3333333333333333 + ], + "shots_by_basis": { + "X": 512, + "Y": 512, + "Z": 512 + } + }, + "jobs": [ + { + "basis": "X", + "task_id": "DF28F46A7519E69A1B1DF8307414D621", + "shots_requested": 512, + "submitted_at_beijing": "2026-07-26T12:31:53+08:00", + "status": "completed", + "originir": "QINIT 1\nCREG 1\nRY q[0],(0.90000000)\nRZ q[0],(0.70000000)\nH q[0]\nMEASURE q[0],c[0]\n", + "counts_before_reconciliation": { + "0": 486, + "1": 26 + }, + "count_reconciliation_delta": 0, + "counts": { + "0": 486, + "1": 26 + }, + "raw_payload": [ + [ + { + "value": [ + 0.9484847, + 0.0515153 + ], + "key": [ + "0x0", + "0x1" + ] + } + ] + ], + "completed_at_beijing": "2026-07-26T12:32:03+08:00" + }, + { + "basis": "Y", + "task_id": "134099D5F971D3B25B49213513CD11BA", + "shots_requested": 512, + "submitted_at_beijing": "2026-07-26T12:32:29+08:00", + "status": "completed", + "originir": "QINIT 1\nCREG 1\nRY q[0],(0.90000000)\nRZ q[0],(0.70000000)\nRX q[0],(1.57079633)\nMEASURE q[0],c[0]\n", + "counts_before_reconciliation": { + "0": 469, + "1": 43 + }, + "count_reconciliation_delta": 0, + "counts": { + "0": 469, + "1": 43 + }, + "raw_payload": [ + [ + { + "value": [ + 0.9164126, + 0.0835874 + ], + "key": [ + "0x0", + "0x1" + ] + } + ] + ], + "completed_at_beijing": "2026-07-26T12:32:36+08:00" + }, + { + "basis": "Z", + "task_id": "D2CEDD69C365796F932599738FC3C0C6", + "shots_requested": 512, + "submitted_at_beijing": "2026-07-26T12:33:00+08:00", + "status": "completed", + "originir": "QINIT 1\nCREG 1\nRY q[0],(0.90000000)\nRZ q[0],(0.70000000)\nMEASURE q[0],c[0]\n", + "counts_before_reconciliation": { + "0": 495, + "1": 17 + }, + "count_reconciliation_delta": 0, + "counts": { + "0": 495, + "1": 17 + }, + "raw_payload": [ + [ + { + "value": [ + 0.9676203, + 0.0323797 + ], + "key": [ + "0x0", + "0x1" + ] + } + ] + ], + "completed_at_beijing": "2026-07-26T12:33:14+08:00" + } + ], + "analysis": { + "analysed_at_beijing": "2026-07-26T12:33:30+08:00", + "family_confidence": 0.95, + "observables": { + "X": { + "hardware": { + "value": 0.8984375, + "standard_error": 0.03779604152268068, + "lower": 0.7157765395198192, + "upper": 1.0810984604801808, + "half_width": 0.18266096048018077, + "shots": 1536, + "empirical_variance": 2.1942385993485343, + "value_range": 6.0, + "failure_probability": 0.01666666666666668, + "minimum_term_matches": 512, + "ideal": 0.5991214669182833, + "absolute_error_from_ideal": 0.2993160330817167, + "ideal_in_statistical_interval": false + }, + "cpuqvm": { + "value": 0.60546875, + "standard_error": 0.04141949998497684, + "lower": 0.4111304033651358, + "upper": 0.7998070966348643, + "half_width": 0.19433834663486424, + "shots": 1536, + "empirical_variance": 2.635123167752443, + "value_range": 6.0, + "failure_probability": 0.01666666666666668, + "minimum_term_matches": 512, + "ideal": 0.5991214669182833, + "absolute_error_from_ideal": 0.006347283081716726, + "ideal_in_statistical_interval": true + } + }, + "Y": { + "hardware": { + "value": 0.83203125, + "standard_error": 0.03877374442445291, + "lower": 0.6462194284164429, + "upper": 1.0178430715835571, + "half_width": 0.18581182158355714, + "shots": 1536, + "empirical_variance": 2.3092274022801305, + "value_range": 6.0, + "failure_probability": 0.01666666666666668, + "minimum_term_matches": 512, + "ideal": 0.5046330500712651, + "absolute_error_from_ideal": 0.32739819992873487, + "ideal_in_statistical_interval": false + }, + "cpuqvm": { + "value": 0.46484375, + "standard_error": 0.04258672154760506, + "lower": 0.2667437769916833, + "upper": 0.6629437230083167, + "half_width": 0.19809997300831672, + "shots": 1536, + "empirical_variance": 2.785733916938111, + "value_range": 6.0, + "failure_probability": 0.01666666666666668, + "minimum_term_matches": 512, + "ideal": 0.5046330500712651, + "absolute_error_from_ideal": 0.03978930007126513, + "ideal_in_statistical_interval": true + } + }, + "Z": { + "hardware": { + "value": 0.93359375, + "standard_error": 0.03723683082657714, + "lower": 0.7527349681036721, + "upper": 1.1144525318963279, + "half_width": 0.18085878189632787, + "shots": 1536, + "empirical_variance": 2.1297892915309444, + "value_range": 6.0, + "failure_probability": 0.01666666666666668, + "minimum_term_matches": 512, + "ideal": 0.6216099682706644, + "absolute_error_from_ideal": 0.3119837817293356, + "ideal_in_statistical_interval": false + }, + "cpuqvm": { + "value": 0.640625, + "standard_error": 0.041073538636694856, + "lower": 0.44740158935806107, + "upper": 0.8338484106419389, + "half_width": 0.19322341064193893, + "shots": 1536, + "empirical_variance": 2.5912866449511403, + "value_range": 6.0, + "failure_probability": 0.01666666666666668, + "minimum_term_matches": 512, + "ideal": 0.6216099682706644, + "absolute_error_from_ideal": 0.019015031729335607, + "ideal_in_statistical_interval": true + } + } + }, + "disclosure": "Empirical Bernstein intervals quantify finite-shot statistical uncertainty only; they do not cover gate noise, readout error, calibration drift, or other systematic bias.", + "readout_mitigation": false, + "hardware_bloch_vector": { + "X": 0.8984375, + "Y": 0.83203125, + "Z": 0.93359375 + }, + "hardware_bloch_norm": 1.5398257149501937, + "validation_outcome": { + "cloud_authentication_and_device_discovery": "passed", + "three_job_budget_and_no_retry_guard": "passed", + "qshadow_ingestion_and_estimation_pipeline": "passed", + "hardware_statistical_accuracy_against_ideal": "not demonstrated" + }, + "interpretation": "The run validates authenticated cloud integration, q0-first count ingestion, per-shot propensity records, simultaneous X/Y/Z estimation with family-wise empirical-Bernstein intervals, and the three-job/1536-shot/no-retry budget. All three hardware intervals exclude their ideal values and the estimated Bloch-vector norm exceeds one, so hardware accuracy is not demonstrated. The run used platform is_amend=True, while platform amendment, gate/readout noise, mapping, and calibration drift are outside QShadow statistical intervals. Because the symmetric X/Y/Z target used a balanced one-third plan, this run does not compare the observable-aware planner against uniform hardware sampling.", + "qshadow_builtin_readout_mitigation": false, + "platform_is_amend": true + }, + "completed_at_beijing": "2026-07-26T12:33:30+08:00", + "runtime_submission_configuration": { + "is_amend": true, + "is_mapping": true, + "is_optimization": true, + "specified_block": null, + "automatic_retry": false + }, + "credential_handling": { + "source": "Windows User environment variable ORIGINQ_API_KEY", + "recorded_in_evidence": false, + "printed_or_logged": false + } +} diff --git a/pyqpanda-algorithm/example/QShadow/qshadow_wukong_validation_20260726.md b/pyqpanda-algorithm/example/QShadow/qshadow_wukong_validation_20260726.md new file mode 100644 index 0000000..1b4fb0c --- /dev/null +++ b/pyqpanda-algorithm/example/QShadow/qshadow_wukong_validation_20260726.md @@ -0,0 +1,56 @@ +# QShadow WK_C180 最小真机验证（2026-07-26） + +本文件记录 QShadow 在本源悟空 `WK_C180` 上的最小端到端验证。三个任务产生的真实 counts 被转换为 q0-first 数据，并按逐 shot 采样概率送入同一个 `ShadowDataset`，随后完成 X/Y/Z 的同时估计、经验 Bernstein 区间和固定任务预算检查。这不是单纯的云接口连通性测试。 + +机器可读原始证据见 [`qshadow_wukong_validation_20260726.json`](qshadow_wukong_validation_20260726.json)。JSON 包含非敏感的 task ID、OriginIR、原始概率、counts、时间戳和统计结果，不包含 API key。 + +## 固定预算与实验设计 + +- 后端：`WK_C180`（Origin Wukong 180） +- SDK：`qpanda3_runtime==1.0.0`、`pyqpanda3==0.3.4` +- 服务：`https://qpanda3-runtime.qpanda.cn` +- 目标态：`RY(0.9) -> RZ(0.7) |0>` +- 可观测量：`X`、`Y`、`Z` +- 计划：X/Y/Z 各 512 shots，逐 shot 采样概率记录为 `(1/3, 1/3, 1/3)` +- 硬限制：3 个任务、总计 1536 shots +- 自动重试：关闭；未切换后端 +- Runtime 配置：`is_amend=True`、`is_mapping=True`、`is_optimization=True` +- QShadow 内建读出误差缓解：无 + +## 云任务证据 + +| 基 | task ID | shots | counts | 计数校准差值 | +|---|---|---:|---|---:| +| X | `DF28F46A7519E69A1B1DF8307414D621` | 512 | `0: 486, 1: 26` | 0 | +| Y | `134099D5F971D3B25B49213513CD11BA` | 512 | `0: 469, 1: 43` | 0 | +| Z | `D2CEDD69C365796F932599738FC3C0C6` | 512 | `0: 495, 1: 17` | 0 | + +三个任务均只提交一次并正常完成。原始概率转换后的 counts 已各自严格求和为 512，无需尾差修正。 + +## 估计结果 + +以下区间为三个可观测量的 95% family-wise 同时经验 Bernstein 区间： + +| 可观测量 | 理想值 | CPUQVM 估计 | WK_C180 估计 | WK_C180 同时区间 | 理想值在区间内 | +|---|---:|---:|---:|---:|---| +| X | 0.599121 | 0.605469 | 0.898438 | [0.715777, 1.081098] | 否 | +| Y | 0.504633 | 0.464844 | 0.832031 | [0.646219, 1.017843] | 否 | +| Z | 0.621610 | 0.640625 | 0.933594 | [0.752735, 1.114453] | 否 | + +由三个独立估计构成的硬件 Bloch 向量范数为 `1.539825715 > 1`，因此不能把这些结果解释成一个物理单量子比特态的无系统偏差重建。 + +## 这次运行验证了什么 + +- **真实数据管线**：完成鉴权、设备发现、OriginIR 转换、任务执行、counts 归一化和 q0-first 摄取，真实 counts 最终进入 QShadow 估计器，而不是停留在云接口层。 +- **同时估计与区间**：同一个含 1536 shots 的 `ShadowDataset` 记录每个 shot 的 `(1/3, 1/3, 1/3)` 倾向概率，并同时产生 X/Y/Z 点估计和 family-wise 置信区间。 +- **任务预算约束**：运行严格保持 3 jobs、1536 shots、0 retry，未切换后端，验证了付费后端场景下的任务数和 shots 边界。 +- **统计边界诊断**：三个硬件区间都排除了理想值，且 Bloch 向量范数大于 1。这说明有限 shots 的区间没有掩盖统计模型之外的设备或平台偏差；它不能定位偏差来源，但能明确提示结果不应按理想态精度解释。 + +## 没有由这次运行验证的结论 + +- 本次 X/Y/Z 目标完全对称，测量计划也是均匀的 `1/3, 1/3, 1/3`，因此没有形成“可观测量感知偏置规划”与均匀规划的真机对照。规划器的代理损失下降和概率下限行为由 CPUQVM 示例及单元测试验证。 +- 三个硬件区间均未覆盖理想值，所以本次运行没有证明 WK_C180 上的估计精度。 +- Runtime 使用平台选项 `is_amend=True`，QShadow 本身未执行读出误差缓解。经验 Bernstein 区间只覆盖有限 shots 的统计不确定性，不覆盖平台修正、门噪声、读出误差、映射或校准漂移。 +- 区间未裁剪到物理范围 `[-1, 1]`，以保留统计公式给出的实际半宽。 + +因此，这三个任务为 QShadow 提供了真实后端上的**数据摄取、同时估计、置信区间和预算控制**验证；硬件精度和偏置规划器的真机优势属于不同的验证问题，不由这组数据作出结论。 diff --git a/pyqpanda-algorithm/pyqpanda_alg/QShadow/README.md b/pyqpanda-algorithm/pyqpanda_alg/QShadow/README.md new file mode 100644 index 0000000..296791a --- /dev/null +++ b/pyqpanda-algorithm/pyqpanda_alg/QShadow/README.md @@ -0,0 +1,272 @@ +# QShadow：可观测量感知的经典阴影估计 + +QShadow 是面向 PyQPanda3 的多可观测量测量与后处理工具。它使用局域 Pauli +`X/Y/Z` 测量，从同一批量子态副本中同时估计多个 Pauli 可观测量；在给定目标 +可观测量后，可优化每个量子比特的测量基概率，并用同时置信区间控制额外 shots。 + +> **准确定位**：经典阴影、局域偏置经典阴影和 Adaptive Pauli Shadows 都是已有 +> 学术成果，本模块不声称发明这些算法。QShadow 的原创贡献范围仅限于本文末尾 +> “原创性与引用边界”中列出的独立工程与算法设计。 + +## 功能 + +- `PauliTerm` / `PauliObservable`：显式 **q0-first** 的 Pauli 数据模型； +- `optimize_basis_probabilities`：以可观测量系数为权重、带严格概率下限的局域测量 + 概率优化； +- `MeasurementPlan` / `ShadowDataset`：记录每个 shot 的测量基、结果和完整采样概率， + 可无偏合并不同测量轮次； +- `estimate_observables`：逆倾向（inverse-propensity）无偏估计与 Bonferroni 同时 + 经验 Bernstein 区间； +- `ShadowBudgetController`：根据目标半宽、最小/最大 shots 和批大小给出下一批预算， + **不会自动提交付费任务**； +- `simulate_pauli_measurements`：不依赖量子后端的 NumPy 小规模参考路径； +- `PyQPandaRunner`：CPUQVM 或用户提供的远程 executor；内置端序转换与 + `max_jobs` 硬上限。 + +## 快速开始 + +```python +import numpy as np +from pyqpanda_alg.QShadow import ( + PauliObservable, + estimate_observables, + optimize_basis_probabilities, + sample_measurement_plan, + simulate_pauli_measurements, +) + +observables = [ + PauliObservable.from_terms( + {"ZZ": -1.0, "XI": 0.25, "IX": 0.25}, name="energy" + ), + PauliObservable.from_terms({"ZI": 1.0}, name="magnetisation_q0"), +] + +optimised = optimize_basis_probabilities( + observables, probability_floor=0.03 +) +plan = sample_measurement_plan(optimised.probabilities, shots=5000, seed=7) + +# 小规模独立参考：|00> +dataset = simulate_pauli_measurements([1, 0, 0, 0], plan, seed=8) +results = estimate_observables(dataset, observables, family_confidence=0.95) +for name, result in results.items(): + print(name, result.value, (result.lower, result.upper)) +``` + +Pauli 字符串的第 `q` 个字符作用于量子比特 `q`，所以 `"XI"` 是 `X(q0)`。 +`ShadowDataset` 内的 bitstring 也采用 q0-first；`PyQPandaRunner` 会将 PyQPanda +通常使用的 MSB-first 计数字符串转换为该约定。 + +## PyQPanda3 本地执行 + +```python +from pyqpanda3.core import CNOT, H, QCircuit +from pyqpanda_alg.QShadow import PyQPandaRunner, run_measurement_plan + + +def prepare(qubits): + return QCircuit() << H(qubits[0]) << CNOT(qubits[0], qubits[1]) + +runner = PyQPandaRunner(prepare, n_qubits=2) # 默认 CPUQVM +dataset = run_measurement_plan(plan, runner, max_jobs=9) +``` + +`max_jobs` 在调用任何 runner **之前**检查，适合为远程后端设置硬任务数上限。 + +## 本源量子云 / 真机接入（后端无关，不保存密钥） + +QShadow 不绑定云 SDK、服务端点或环境变量名称，也不接收和保存 API 密钥。应用层应使用 +当前受支持的 SDK 在 QShadow 之外完成鉴权，并把一个 +``executor(program, shots) -> counts`` 回调交给 ``PyQPandaRunner``。凭据只应从本地 +环境变量或密钥存储读取，不得写入源码、Notebook、PR、日志或命令行参数。 + +下面是**接口合约伪代码**，其中 ``submit_once_and_wait`` 由应用层适配器实现；它不是 +某个云 SDK 的可直接运行示例，也不应隐式重试或切换后端： + +```python +from pyqpanda_alg.QShadow import PyQPandaRunner, run_measurement_plan + +job_ids = [] + + +def execute(program, shots): + result = submit_once_and_wait(program, shots=shots) # 应用层 SDK 适配器 + job_ids.append(result.task_id) # 任务 ID 非凭据 + return result.counts # Mapping[str, int] + +required_jobs = len(plan.grouped_shots) +if required_jobs > 3: + raise RuntimeError("remote job budget exceeded") + +runner = PyQPandaRunner(prepare, 2, executor=execute) +dataset = run_measurement_plan(plan, runner, max_jobs=3) +print(job_ids) +``` + +``max_jobs`` 会在任何 executor 调用前检查压缩后的唯一测量基数量。已验证的可选 +``qpanda3_runtime==1.0.0`` 路径及非敏感原始证据见下文；该 SDK 不是 QShadow 的 +核心依赖，服务地址和鉴权方式应以部署时的官方文档为准。 + +建议的真机证据链： + +1. 固定代码版本、随机种子、后端名称、芯片/映射信息和 UTC 时间； +2. 先通过 NumPy 与 CPUQVM，再提交 1–3 比特、受控任务数和 shots 的真机实验； +3. 保存每个任务 ID、原始 counts、测量计划和采样概率； +4. 报告估计值、统计区间及理想模拟器偏差，不把统计区间误称为硬件噪声区间； +5. 如未做读出误差缓解，应明确披露。QShadow 当前不会隐式修正设备噪声。 + +### 已执行的 WK_C180 最小验证（2026-07-26） + +在明确上限后，验证通过 ``qpanda3_runtime==1.0.0`` 向 ``WK_C180`` 提交了 +X/Y/Z 各 512 shots，合计恰好 3 个任务、1536 shots，且没有自动重试。非敏感的 +task ID、OriginIR、原始概率、counts 和分析保存在 +[`qshadow_wukong_validation_20260726.json`](../../example/QShadow/qshadow_wukong_validation_20260726.json)， +解释说明见 +[`qshadow_wukong_validation_20260726.md`](../../example/QShadow/qshadow_wukong_validation_20260726.md)。 + +这三个任务不只是云接口连通性测试：WK_C180 的原始 counts 经归一化和 q0-first +转换后进入同一个 ``ShadowDataset``，逐 shot 记录 ``(1/3, 1/3, 1/3)`` 概率， +并产生 X/Y/Z 的同时点估计与经验 Bernstein 区间；运行同时保持了 3 jobs、 +1536 shots 和 0 retry 的硬预算。 + +硬件 X/Y/Z 估计 ``0.898438 / 0.832031 / 0.933594`` 的同时区间均未覆盖理想值 +``0.599121 / 0.504633 / 0.621610``，独立 Bloch 向量范数为 +``1.539826 > 1``。这使统计模型之外的设备或平台偏差在结果中可见，而没有被有限 +shots 的区间掩盖。本次 Runtime 使用 ``is_amend=True``，QShadow 本身没有内建 +读出缓解，因此不能据此声称硬件精度成功。 + +本次目标集合 X/Y/Z 完全对称，使用的是均匀计划，未比较偏置规划器与均匀规划的 +真机表现。规划器的代理损失和概率下限行为由 CPUQVM 示例与测试验证；本记录验证的 +是真实 counts 摄取、同时估计、置信区间和任务预算流程。 + +## 方法 + +### 无偏局域 Pauli 阴影估计 + +对 Pauli 项 `P`，第 `t` 次测量在量子比特 `q` 选择基 +`B[t,q] ~ beta[t,q]`，得到本征值 `s[t,q] in {-1,+1}`。QShadow 使用 + +```text +X_t(P) = 1[B_t covers P] * product(s[t,q]) + / product(beta[t,q,P[q]]) +``` + +作为 `

` 的无偏估计。分母保留每个 shot 当时的完整概率，因此先均匀测量、再 +使用偏置测量的多轮数据可以直接合并。 + +### 可观测量感知规划 + +对目标 Pauli 项 `P_j`、系数 `c_j` 和用户权重 `w_j`，优化器最小化对角二阶矩代理 + +```text +L(beta) = sum_j w_j c_j^2 + / product(beta[q, P_j[q]], q in support(P_j)). +``` + +固定其余量子比特后，单个量子比特的子问题为 +`min sum_a A_a / beta_a`；无下限时解析解满足 `beta_a ∝ sqrt(A_a)`。本实现加入 +`beta_a >= probability_floor`，用活动集求出每个块的精确最优点并循环至收敛。 +该目标是**状态无关的对角代理**，不是完整协方差，也不保证最小化未知量子态上的 +真实方差。 + +### 可复现离线对照 + +[`qshadow_planning_benchmark.py`](../../example/QShadow/qshadow_planning_benchmark.py) +在相同三比特状态、目标可观测量和 shots 下比较均匀规划与 QShadow 规划。以下结果 +使用固定种子、1000 轮 × 2000 shots；顺序预算部分另运行 100 轮： + +```bash +PYTHONPATH=pyqpanda-algorithm python \ + pyqpanda-algorithm/example/QShadow/qshadow_planning_benchmark.py \ + --trials 1000 --shots 2000 --budget-trials 100 \ + --target-half-width 0.22 --seed 2000000 +``` + +- 规划代理损失：``24.54 -> 5.964549``； +- 聚合 RMSE 下降 ``60.20%``；两个可观测量的 RMSE 分别下降 ``61.72%`` 和 + ``57.68%``； +- 平均同时区间半宽分别下降 ``62.35%`` 和 ``53.42%``； +- 达到相同 ``0.22`` 目标半宽时，平均总 shots 从 ``7284.9`` 降至 ``1892.8``， + 减少 ``74.02%``，且 100 轮均无硬上限违规； +- 本实验两方案的单项和 family-wise 经验覆盖率均为 ``1.000``，说明区间在该配置下 + 较保守；这不是“精确校准为 95%”的声明。 + +1000 轮中优化方案两个偏差的 z-score 为 ``0.28`` 和 ``-0.26``，与零一致。上述 +改善只针对固定状态、可观测量、概率下限和种子日程，不推广为所有问题上的普遍优势。 + +### 置信预算 + +对每个可观测量的逐 shot 样本，QShadow 使用 + +```text +radius = sqrt(2 s^2 log(3/delta) / n) + 3 R log(3/delta) / n +``` + +作为保守经验 Bernstein 半宽，并以 Bonferroni 分配 `delta`，形成多个可观测量的 +同时区间。`R` 来自概率下限下的逆倾向确定性范围。高权重 Pauli 项可能令该范围很 +大，因此区间常比高斯误差条更保守。 + +## 限制 + +- 目前仅支持实系数 Hermitian Pauli 和局域 `X/Y/Z` 测量； +- 不包含读出误差缓解、零噪声外推或完整状态层析； +- 经验 Bernstein 区间控制抽样随机性，不覆盖校准漂移、门噪声等系统偏差； +- 状态向量参考采样是指数复杂度，只用于小规模校验； +- 每个唯一 basis 通常对应一个后端任务，真机前必须检查 + `len(plan.grouped_shots)` 并设置 `max_jobs`； +- “规划损失下降”不等同于所有量子态、所有可观测量上的实际均方误差必然下降。 + +## 原创性与引用边界 + +### 已有工作（明确归因，绝不声称首创） + +1. H.-Y. Huang, R. Kueng, J. Preskill, *Predicting many properties of a + quantum system from very few measurements*, Nature Physics 16, 1050–1057 + (2020), [arXiv:2002.08953](https://arxiv.org/abs/2002.08953)。 +2. C. Hadfield, S. Bravyi, R. Raymond, A. Mezzacapo, *Measurements of Quantum + Hamiltonians with Locally-Biased Classical Shadows*, + [arXiv:2006.15788](https://arxiv.org/abs/2006.15788)。 +3. C. Hadfield, *Adaptive Pauli Shadows for Energy Estimation*, + [arXiv:2105.12207](https://arxiv.org/abs/2105.12207)。这是已有自适应方法； + QShadow 不把“自适应 Pauli 阴影”作为自己的发明。 + +### 本贡献可主张的原创范围 + +- 面向 PyQPanda3、q0-first 且后端无关的完整数据与 runner API； +- 基于上述明确代理目标、带概率下限的活动集块坐标实现及收敛诊断； +- 支持异构测量轮次的逐 shot 概率记录与无偏合并； +- 同时经验 Bernstein 区间到硬 shots 上限之间的显式预算控制； +- `max_jobs` 付费后端安全阈值、端序回归测试、NumPy/CPUQVM 双参考路径； +- 本仓库中的代码、测试、示例和文档均独立编写，未复制论文作者或其他框架的 + 源文件。 + +## 测试 + +```bash +PYTHONPATH=pyqpanda-algorithm python -m pytest test/QShadow -q -o addopts="" +``` + +测试覆盖 Pauli 代数、规划损失单调性、概率下限、多轮无偏合并、置信预算、计数 +校验、任务上限、端序，以及 CPUQVM 的 X/Y/Z 基符号。 + +--- + +## English summary + +QShadow is a PyQPanda3-oriented workflow for estimating many Pauli observables +from local X/Y/Z measurements. It provides observable-aware locally biased +planning, per-shot propensity records, unbiased post-processing, simultaneous +empirical-Bernstein intervals, explicit shot-budget decisions, and local/cloud +runner adapters. Classical shadows and locally/adaptively biased Pauli shadows +are prior work cited above; this contribution claims only its independently +written PyQPanda integration, constrained coordinate optimiser, data model, +budget workflow, examples, documentation, and tests. Device noise is not +silently mitigated, and real-hardware results must report backend metadata, job +IDs, raw counts, statistical uncertainty, and simulator discrepancy. A disclosed +WK_C180 run on 2026-07-26 used exactly three 512-shot jobs with no retry. Cloud +integration and budget controls passed, while hardware accuracy was not +demonstrated: all three simultaneous intervals missed the ideal values and the +independently estimated Bloch-vector norm exceeded one. The run used the +platform option ``is_amend=True``; full raw evidence and limitations are stored +under ``example/QShadow/qshadow_wukong_validation_20260726.*``. diff --git a/pyqpanda-algorithm/pyqpanda_alg/QShadow/__init__.py b/pyqpanda-algorithm/pyqpanda_alg/QShadow/__init__.py new file mode 100644 index 0000000..d179d4e --- /dev/null +++ b/pyqpanda-algorithm/pyqpanda_alg/QShadow/__init__.py @@ -0,0 +1,61 @@ +"""QShadow: observable-aware local-Pauli classical-shadow estimation. + +QShadow combines independently implemented locally-biased measurement planning, +backend-neutral records, unbiased inverse-propensity estimation, simultaneous +empirical-Bernstein intervals, and explicit shot-budget control. See README.md +for references, originality scope, and limitations. +""" + +from .dataset import ShadowDataset +from .estimator import ( + BudgetReport, + ShadowBudgetController, + ShadowEstimate, + empirical_bernstein_radius, + estimate_observable, + estimate_observables, + observable_shadow_samples, + term_shadow_samples, +) +from .observable import PauliObservable, PauliTerm, pauli_expectation +from .planner import ( + AXES, + BasisOptimizationResult, + MeasurementPlan, + optimize_basis_probabilities, + planning_loss, + sample_measurement_plan, + uniform_probabilities, +) +from .sampler import ( + PyQPandaRunner, + build_pyqpanda_program, + run_measurement_plan, + simulate_pauli_measurements, +) + +__all__ = [ + "AXES", + "BasisOptimizationResult", + "BudgetReport", + "MeasurementPlan", + "PauliObservable", + "PauliTerm", + "PyQPandaRunner", + "ShadowBudgetController", + "ShadowDataset", + "ShadowEstimate", + "build_pyqpanda_program", + "empirical_bernstein_radius", + "estimate_observable", + "estimate_observables", + "observable_shadow_samples", + "optimize_basis_probabilities", + "pauli_expectation", + "planning_loss", + "run_measurement_plan", + "sample_measurement_plan", + "simulate_pauli_measurements", + "term_shadow_samples", + "uniform_probabilities", +] diff --git a/pyqpanda-algorithm/pyqpanda_alg/QShadow/dataset.py b/pyqpanda-algorithm/pyqpanda_alg/QShadow/dataset.py new file mode 100644 index 0000000..3740838 --- /dev/null +++ b/pyqpanda-algorithm/pyqpanda_alg/QShadow/dataset.py @@ -0,0 +1,148 @@ +"""Backend-neutral measurement records for local-Pauli classical shadows.""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Sequence + +import numpy as np + +from .planner import AXES, _validate_probabilities + +_AXIS_TO_INDEX = {axis: index for index, axis in enumerate(AXES)} + + +@dataclass(frozen=True, slots=True) +class ShadowDataset: + """Local-Pauli bases, eigenvalue outcomes, and sampling probabilities. + + Arrays use the shape ``(shots, n_qubits)`` except ``probabilities``, whose + shape is ``(shots, n_qubits, 3)`` in X/Y/Z order. Keeping the complete + distribution for every shot allows datasets from different planning rounds + to be concatenated without bias. + """ + + bases: np.ndarray + outcomes: np.ndarray + probabilities: np.ndarray + + def __post_init__(self) -> None: + bases = np.asarray(self.bases, dtype=np.int8).copy() + outcomes = np.asarray(self.outcomes, dtype=np.int8).copy() + probabilities = np.asarray(self.probabilities, dtype=float).copy() + if bases.ndim != 2 or bases.shape[0] == 0 or bases.shape[1] == 0: + raise ValueError("bases must have shape (positive shots, positive qubits)") + if outcomes.shape != bases.shape: + raise ValueError("outcomes must have the same shape as bases") + if probabilities.shape != (*bases.shape, 3): + raise ValueError( + "probabilities must have shape (shots, n_qubits, 3)" + ) + if np.any((bases < 0) | (bases > 2)): + raise ValueError("basis indices must be 0 (X), 1 (Y), or 2 (Z)") + if np.any((outcomes != 1) & (outcomes != -1)): + raise ValueError("outcomes must be Pauli eigenvalues +1 or -1") + if not np.all(np.isfinite(probabilities)) or np.any(probabilities <= 0): + raise ValueError("all basis probabilities must be finite and positive") + if not np.allclose(probabilities.sum(axis=2), 1.0, atol=1e-10): + raise ValueError("each per-shot X/Y/Z distribution must sum to one") + bases.setflags(write=False) + outcomes.setflags(write=False) + probabilities.setflags(write=False) + object.__setattr__(self, "bases", bases) + object.__setattr__(self, "outcomes", outcomes) + object.__setattr__(self, "probabilities", probabilities) + + @classmethod + def from_strings( + cls, + bases: Sequence[str], + bitstrings: Sequence[str], + probabilities: Sequence[Sequence[float]] | np.ndarray, + ) -> "ShadowDataset": + """Construct records from q0-first basis and outcome strings. + + A bit ``0`` maps to eigenvalue ``+1`` and bit ``1`` to ``-1``. + ``probabilities`` may be one static ``(n_qubits, 3)`` distribution or a + per-shot ``(shots, n_qubits, 3)`` array. + """ + + basis_strings = tuple(bases) + outcome_strings = tuple(bitstrings) + if not basis_strings or len(basis_strings) != len(outcome_strings): + raise ValueError("bases and bitstrings must have equal positive length") + n_qubits = len(basis_strings[0]) + basis_array: np.ndarray = np.empty((len(basis_strings), n_qubits), dtype=np.int8) + outcome_array: np.ndarray = np.empty_like(basis_array) + for shot, (basis, bits) in enumerate(zip(basis_strings, outcome_strings)): + if len(basis) != n_qubits or set(basis) - set(AXES): + raise ValueError("bases must be equal-length q0-first X/Y/Z strings") + if len(bits) != n_qubits or set(bits) - {"0", "1"}: + raise ValueError("bitstrings must be equal-length q0-first binary strings") + basis_array[shot] = [_AXIS_TO_INDEX[axis] for axis in basis] + outcome_array[shot] = [1 if bit == "0" else -1 for bit in bits] + + raw_probabilities: np.ndarray = np.asarray(probabilities, dtype=float) + if raw_probabilities.ndim == 2: + static = _validate_probabilities( + raw_probabilities, n_qubits=n_qubits + ) + probability_array = np.broadcast_to( + static, (len(basis_strings), n_qubits, 3) + ).copy() + elif raw_probabilities.shape == (len(basis_strings), n_qubits, 3): + probability_array = raw_probabilities.copy() + else: + raise ValueError( + "probabilities must have shape (n_qubits, 3) or " + "(shots, n_qubits, 3)" + ) + return cls(basis_array, outcome_array, probability_array) + + @classmethod + def concatenate(cls, *datasets: "ShadowDataset") -> "ShadowDataset": + """Combine independently collected rounds with compatible qubit counts.""" + + if not datasets: + raise ValueError("at least one dataset is required") + if not all(isinstance(dataset, cls) for dataset in datasets): + raise TypeError("all inputs must be ShadowDataset instances") + n_qubits = datasets[0].n_qubits + if any(dataset.n_qubits != n_qubits for dataset in datasets): + raise ValueError("all datasets must use the same number of qubits") + return cls( + np.concatenate([dataset.bases for dataset in datasets], axis=0), + np.concatenate([dataset.outcomes for dataset in datasets], axis=0), + np.concatenate([dataset.probabilities for dataset in datasets], axis=0), + ) + + @property + def shots(self) -> int: + return int(self.bases.shape[0]) + + @property + def n_qubits(self) -> int: + return int(self.bases.shape[1]) + + @property + def basis_strings(self) -> tuple[str, ...]: + return tuple("".join(AXES[index] for index in row) for row in self.bases) + + @property + def bitstrings(self) -> tuple[str, ...]: + return tuple( + "".join("0" if value == 1 else "1" for value in row) + for row in self.outcomes + ) + + def take(self, indices: Sequence[int] | slice) -> "ShadowDataset": + """Return a validated shot subset.""" + + bases = self.bases[indices] + outcomes = self.outcomes[indices] + probabilities = self.probabilities[indices] + if bases.ndim == 1: + bases = bases[None, :] + outcomes = outcomes[None, :] + probabilities = probabilities[None, :, :] + return ShadowDataset(bases, outcomes, probabilities) diff --git a/pyqpanda-algorithm/pyqpanda_alg/QShadow/estimator.py b/pyqpanda-algorithm/pyqpanda_alg/QShadow/estimator.py new file mode 100644 index 0000000..b0eab25 --- /dev/null +++ b/pyqpanda-algorithm/pyqpanda_alg/QShadow/estimator.py @@ -0,0 +1,325 @@ +"""Unbiased Pauli-shadow estimators and confidence-budget control.""" + +from __future__ import annotations + +from dataclasses import dataclass +from math import ceil, log, sqrt +from typing import Mapping, Sequence + +import numpy as np + +from .dataset import ShadowDataset +from .observable import PauliObservable, PauliTerm, coerce_observables +from .planner import _AXIS_TO_INDEX + + +def term_shadow_samples(dataset: ShadowDataset, term: PauliTerm) -> np.ndarray: + """Return one inverse-propensity shadow estimate per shot for ``term``.""" + + if term.n_qubits != dataset.n_qubits: + raise ValueError("term and dataset qubit counts differ") + if term.is_identity: + return np.full(dataset.shots, term.coefficient, dtype=float) + + support: np.ndarray = np.asarray(term.support, dtype=int) + required_axes: np.ndarray = np.asarray([_AXIS_TO_INDEX[term.pauli[qubit]] for qubit in support], dtype=int) + matches: np.ndarray = np.all(dataset.bases[:, support] == required_axes, axis=1) + coverage: np.ndarray = np.ones(dataset.shots, dtype=float) + signs: np.ndarray = np.ones(dataset.shots, dtype=float) + for qubit, axis in zip(support, required_axes): + coverage *= dataset.probabilities[:, qubit, axis] + signs *= dataset.outcomes[:, qubit] + samples: np.ndarray = np.zeros(dataset.shots, dtype=float) + samples[matches] = term.coefficient * signs[matches] / coverage[matches] + return samples + + +def observable_shadow_samples( + dataset: ShadowDataset, observable: PauliObservable +) -> np.ndarray: + """Return per-shot unbiased estimates of a Pauli observable.""" + + if observable.n_qubits != dataset.n_qubits: + raise ValueError("observable and dataset qubit counts differ") + samples: np.ndarray = np.zeros(dataset.shots, dtype=float) + for term in observable.terms: + samples += term_shadow_samples(dataset, term) + return samples + + +def _observable_value_range( + dataset: ShadowDataset, observable: PauliObservable +) -> float: + """Conservative common range width for all per-shot estimators.""" + + absolute_bound: np.ndarray = np.zeros(dataset.shots, dtype=float) + for term in observable.terms: + if term.is_identity or term.coefficient == 0.0: + continue + inverse_coverage: np.ndarray = np.ones(dataset.shots, dtype=float) + for qubit in term.support: + axis = _AXIS_TO_INDEX[term.pauli[qubit]] + inverse_coverage /= dataset.probabilities[:, qubit, axis] + absolute_bound += abs(term.coefficient) * inverse_coverage + return float(2.0 * np.max(absolute_bound, initial=0.0)) + + +def empirical_bernstein_radius( + samples: Sequence[float] | np.ndarray, + *, + value_range: float, + failure_probability: float, +) -> float: + """Return a conservative two-sided empirical Bernstein half-width. + + The implemented form is + + ``sqrt(2 s^2 log(3/delta) / n) + 3 R log(3/delta) / n``, + + where ``s^2`` is the unbiased sample variance and ``R`` is a known range + width. The result is capped by ``R``, the deterministic worst-case mean + deviation. It is intentionally conservative for high-weight Pauli sums. + """ + + values: np.ndarray = np.asarray(samples, dtype=float) + if values.ndim != 1 or values.size == 0: + raise ValueError("samples must be a non-empty one-dimensional sequence") + if not np.all(np.isfinite(values)): + raise ValueError("samples must be finite") + if not np.isfinite(value_range) or value_range < 0: + raise ValueError("value_range must be finite and non-negative") + if not 0.0 < failure_probability < 1.0: + raise ValueError("failure_probability must lie between zero and one") + if value_range == 0.0: + return 0.0 + if values.size == 1: + return float(value_range) + variance = float(np.var(values, ddof=1)) + logarithm = log(3.0 / failure_probability) + radius = sqrt(2.0 * variance * logarithm / values.size) + radius += 3.0 * value_range * logarithm / values.size + return float(min(value_range, radius)) + + +@dataclass(frozen=True, slots=True) +class ShadowEstimate: + """Point estimate, diagnostics, and a finite-sample confidence interval.""" + + observable: str + value: float + standard_error: float + lower: float + upper: float + shots: int + empirical_variance: float + value_range: float + failure_probability: float + minimum_term_matches: int + + @property + def half_width(self) -> float: + return (self.upper - self.lower) / 2.0 + + @property + def confidence(self) -> float: + return 1.0 - self.failure_probability + + +def _minimum_term_matches( + dataset: ShadowDataset, observable: PauliObservable +) -> int: + counts: list[int] = [] + for term in observable.terms: + if term.is_identity or term.coefficient == 0.0: + continue + matched: np.ndarray = np.ones(dataset.shots, dtype=bool) + for qubit in term.support: + matched &= dataset.bases[:, qubit] == _AXIS_TO_INDEX[term.pauli[qubit]] + counts.append(int(matched.sum())) + return min(counts, default=dataset.shots) + + +def estimate_observable( + dataset: ShadowDataset, + observable: PauliObservable, + *, + confidence: float = 0.95, + failure_probability: float | None = None, +) -> ShadowEstimate: + """Estimate one observable from any compatible QShadow dataset.""" + + if not 0.0 < confidence < 1.0: + raise ValueError("confidence must lie between zero and one") + delta = 1.0 - confidence if failure_probability is None else failure_probability + if not 0.0 < delta < 1.0: + raise ValueError("failure_probability must lie between zero and one") + samples = observable_shadow_samples(dataset, observable) + value = float(np.mean(samples)) + variance = float(np.var(samples, ddof=1)) if dataset.shots > 1 else 0.0 + standard_error = sqrt(variance / dataset.shots) + value_range = _observable_value_range(dataset, observable) + radius = empirical_bernstein_radius( + samples, value_range=value_range, failure_probability=delta + ) + return ShadowEstimate( + observable=observable.name, + value=value, + standard_error=standard_error, + lower=value - radius, + upper=value + radius, + shots=dataset.shots, + empirical_variance=variance, + value_range=value_range, + failure_probability=delta, + minimum_term_matches=_minimum_term_matches(dataset, observable), + ) + + +def estimate_observables( + dataset: ShadowDataset, + observables: PauliObservable | Sequence[PauliObservable], + *, + family_confidence: float = 0.95, +) -> dict[str, ShadowEstimate]: + """Estimate many observables with Bonferroni simultaneous confidence.""" + + checked = coerce_observables(observables) + if any(observable.n_qubits != dataset.n_qubits for observable in checked): + raise ValueError("observable and dataset qubit counts differ") + if not 0.0 < family_confidence < 1.0: + raise ValueError("family_confidence must lie between zero and one") + per_observable_delta = (1.0 - family_confidence) / len(checked) + return { + observable.name: estimate_observable( + dataset, + observable, + failure_probability=per_observable_delta, + ) + for observable in checked + } + + +def _required_total_shots( + *, + variance: float, + value_range: float, + failure_probability: float, + target_half_width: float, +) -> int: + if value_range == 0.0: + return 1 + logarithm = log(3.0 / failure_probability) + linear = sqrt(max(0.0, 2.0 * variance * logarithm)) + constant = 3.0 * value_range * logarithm + root_n = ( + linear + sqrt(linear * linear + 4.0 * target_half_width * constant) + ) / (2.0 * target_half_width) + return max(2, ceil(root_n * root_n)) + + +@dataclass(frozen=True, slots=True) +class BudgetReport: + """Decision returned by :class:`ShadowBudgetController`.""" + + estimates: Mapping[str, ShadowEstimate] + complete: bool + exhausted: bool + next_batch_shots: int + reason: str + current_shots: int + max_shots: int + + +class ShadowBudgetController: + """Translate simultaneous confidence widths into a bounded shot budget. + + The controller does not silently submit jobs. Callers evaluate a dataset, + inspect :class:`BudgetReport`, explicitly collect ``next_batch_shots`` more + records, concatenate them, and evaluate again. This makes paid QPU use + auditable and keeps a hard ``max_shots`` ceiling. + """ + + def __init__( + self, + observables: PauliObservable | Sequence[PauliObservable], + *, + target_half_width: float | Mapping[str, float], + family_confidence: float = 0.95, + min_shots: int = 100, + max_shots: int = 10_000, + batch_shots: int = 200, + ) -> None: + self.observables = coerce_observables(observables) + if isinstance(target_half_width, Mapping): + targets = {name: float(value) for name, value in target_half_width.items()} + expected_names = {observable.name for observable in self.observables} + if set(targets) != expected_names: + raise ValueError("target_half_width mapping must match observable names") + else: + value = float(target_half_width) + targets = {observable.name: value for observable in self.observables} + if any(not np.isfinite(value) or value <= 0 for value in targets.values()): + raise ValueError("target half-widths must be finite and positive") + if not 0.0 < family_confidence < 1.0: + raise ValueError("family_confidence must lie between zero and one") + for name, value in ( + ("min_shots", min_shots), + ("max_shots", max_shots), + ("batch_shots", batch_shots), + ): + if not isinstance(value, int) or isinstance(value, bool) or value <= 0: + raise ValueError(f"{name} must be a positive integer") + if min_shots > max_shots: + raise ValueError("min_shots cannot exceed max_shots") + self.targets = targets + self.family_confidence = family_confidence + self.min_shots = min_shots + self.max_shots = max_shots + self.batch_shots = batch_shots + + def evaluate(self, dataset: ShadowDataset) -> BudgetReport: + """Assess current confidence widths without performing new measurements.""" + + if dataset.n_qubits != self.observables[0].n_qubits: + raise ValueError("dataset and controller qubit counts differ") + estimates = estimate_observables( + dataset, self.observables, family_confidence=self.family_confidence + ) + widths_reached = all( + estimate.half_width <= self.targets[name] + for name, estimate in estimates.items() + ) + complete = dataset.shots >= self.min_shots and widths_reached + exhausted = dataset.shots >= self.max_shots and not complete + if complete: + return BudgetReport( + estimates, True, False, 0, "target_reached", dataset.shots, self.max_shots + ) + if exhausted: + return BudgetReport( + estimates, False, True, 0, "max_shots_reached", dataset.shots, self.max_shots + ) + + required = self.min_shots + for name, estimate in estimates.items(): + required = max( + required, + _required_total_shots( + variance=estimate.empirical_variance, + value_range=estimate.value_range, + failure_probability=estimate.failure_probability, + target_half_width=self.targets[name], + ), + ) + remaining = self.max_shots - dataset.shots + suggested = max(self.batch_shots, required - dataset.shots) + next_batch = min(remaining, suggested) + return BudgetReport( + estimates, + False, + False, + next_batch, + "more_shots_required", + dataset.shots, + self.max_shots, + ) diff --git a/pyqpanda-algorithm/pyqpanda_alg/QShadow/observable.py b/pyqpanda-algorithm/pyqpanda_alg/QShadow/observable.py new file mode 100644 index 0000000..955ad8b --- /dev/null +++ b/pyqpanda-algorithm/pyqpanda_alg/QShadow/observable.py @@ -0,0 +1,226 @@ +"""Pauli observable data structures used by :mod:`pyqpanda_alg.QShadow`. + +The public convention in this module is deliberately explicit: character ``q`` +in a Pauli string acts on qubit ``q``. Therefore, ``"XI"`` means X on qubit 0 +and identity on qubit 1. This avoids leaking backend-specific bit-string +endianness into the statistical code. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Iterable, Mapping, Sequence + +import numpy as np + +_ALLOWED_PAULIS = frozenset("IXYZ") + + +def _normalise_pauli(pauli: str) -> str: + if not isinstance(pauli, str): + raise TypeError("pauli must be a string") + normalised = "".join(pauli.split()).upper() + if not normalised: + raise ValueError("a Pauli string cannot be empty") + invalid = sorted(set(normalised) - _ALLOWED_PAULIS) + if invalid: + raise ValueError( + f"invalid Pauli character(s) {invalid}; expected only I, X, Y, Z" + ) + return normalised + + +@dataclass(frozen=True, slots=True) +class PauliTerm: + """A real-weighted Pauli string. + + Parameters + ---------- + pauli: + Pauli string in q0-first order. Character ``pauli[q]`` acts on qubit + ``q``. + coefficient: + Finite real coefficient. QShadow estimates Hermitian observables, so a + complex coefficient is intentionally rejected. + """ + + pauli: str + coefficient: float = 1.0 + + def __post_init__(self) -> None: + object.__setattr__(self, "pauli", _normalise_pauli(self.pauli)) + if isinstance(self.coefficient, complex): + if self.coefficient.imag != 0: + raise ValueError("PauliTerm coefficient must be real") + coefficient = float(self.coefficient.real) + else: + try: + coefficient = float(self.coefficient) + except (TypeError, ValueError) as exc: + raise TypeError("PauliTerm coefficient must be a real number") from exc + if not np.isfinite(coefficient): + raise ValueError("PauliTerm coefficient must be finite") + object.__setattr__(self, "coefficient", coefficient) + + @property + def n_qubits(self) -> int: + """Number of qubits addressed by the Pauli string.""" + + return len(self.pauli) + + @property + def support(self) -> tuple[int, ...]: + """Qubit indices on which the term is non-identity.""" + + return tuple(q for q, axis in enumerate(self.pauli) if axis != "I") + + @property + def weight(self) -> int: + """Pauli weight (number of non-identity factors).""" + + return len(self.support) + + @property + def is_identity(self) -> bool: + """Whether this term is the all-identity operator.""" + + return self.weight == 0 + + +@dataclass(frozen=True, slots=True) +class PauliObservable: + """A named real linear combination of Pauli strings.""" + + terms: tuple[PauliTerm, ...] + name: str = "observable" + + def __post_init__(self) -> None: + terms = tuple(self.terms) + if not terms: + raise ValueError("an observable must contain at least one Pauli term") + if not all(isinstance(term, PauliTerm) for term in terms): + raise TypeError("terms must contain only PauliTerm instances") + n_qubits = terms[0].n_qubits + if any(term.n_qubits != n_qubits for term in terms): + raise ValueError("all Pauli terms in an observable must have equal length") + if not isinstance(self.name, str) or not self.name.strip(): + raise ValueError("observable name must be a non-empty string") + object.__setattr__(self, "terms", terms) + object.__setattr__(self, "name", self.name.strip()) + + @classmethod + def from_terms( + cls, + terms: ( + Mapping[str, float] + | Iterable[PauliTerm | tuple[str, float]] + ), + *, + name: str = "observable", + ) -> "PauliObservable": + """Build an observable and combine duplicate Pauli strings. + + ``terms`` may be a mapping such as ``{"ZI": 0.5, "IZ": 0.5}``, an + iterable of ``(pauli, coefficient)`` pairs, or existing + :class:`PauliTerm` objects. Duplicate strings are summed in insertion + order. Zero coefficients are retained so the number of qubits remains + unambiguous even for a zero observable. + """ + + source: Iterable[PauliTerm | tuple[str, float]] + source = terms.items() if isinstance(terms, Mapping) else terms + combined: dict[str, float] = {} + for item in source: + term = item if isinstance(item, PauliTerm) else PauliTerm(*item) + combined[term.pauli] = combined.get(term.pauli, 0.0) + term.coefficient + if not combined: + raise ValueError("terms cannot be empty") + return cls( + tuple(PauliTerm(pauli, coefficient) for pauli, coefficient in combined.items()), + name=name, + ) + + @property + def n_qubits(self) -> int: + """Number of qubits addressed by this observable.""" + + return self.terms[0].n_qubits + + def exact_expectation(self, statevector: Sequence[complex]) -> float: + """Return the exact statevector expectation for validation. + + The statevector follows the usual little-endian simulator layout: bit + ``q`` of an amplitude index is qubit ``q``. This method is a reference + calculation for examples and tests, not a replacement for measurements. + """ + + value = sum( + term.coefficient * pauli_expectation(statevector, term.pauli) + for term in self.terms + ) + return float(value) + + +def pauli_expectation(statevector: Sequence[complex], pauli: str) -> float: + """Compute ```` without constructing a dense Pauli matrix.""" + + pauli = _normalise_pauli(pauli) + state = np.asarray(statevector, dtype=np.complex128) + if state.ndim != 1: + raise ValueError("statevector must be one-dimensional") + expected_size = 1 << len(pauli) + if state.size != expected_size: + raise ValueError( + f"statevector has length {state.size}, expected {expected_size} " + f"for {len(pauli)} qubits" + ) + norm = float(np.vdot(state, state).real) + if not np.isfinite(norm) or not np.isclose(norm, 1.0, atol=1e-10): + raise ValueError("statevector must be finite and normalised") + + transformed = np.zeros_like(state) + for source, amplitude in enumerate(state): + target = source + phase = 1.0 + 0.0j + for qubit, axis in enumerate(pauli): + if axis == "I": + continue + bit = (source >> qubit) & 1 + if axis == "X": + target ^= 1 << qubit + elif axis == "Y": + target ^= 1 << qubit + phase *= 1j if bit == 0 else -1j + else: # Z + phase *= 1.0 if bit == 0 else -1.0 + transformed[target] += phase * amplitude + + expectation = np.vdot(state, transformed) + if abs(expectation.imag) > 1e-9: + raise ArithmeticError( + "a Hermitian Pauli expectation acquired a non-negligible imaginary part" + ) + return float(expectation.real) + + +def coerce_observables( + observables: PauliObservable | Sequence[PauliObservable], +) -> tuple[PauliObservable, ...]: + """Validate and normalise one or more same-size observables.""" + + result: tuple[PauliObservable, ...] + if isinstance(observables, PauliObservable): + result = (observables,) + else: + result = tuple(observables) + if not result: + raise ValueError("at least one observable is required") + if not all(isinstance(observable, PauliObservable) for observable in result): + raise TypeError("observables must contain only PauliObservable instances") + n_qubits = result[0].n_qubits + if any(observable.n_qubits != n_qubits for observable in result): + raise ValueError("all observables must act on the same number of qubits") + names = [observable.name for observable in result] + if len(set(names)) != len(names): + raise ValueError("observable names must be unique") + return result diff --git a/pyqpanda-algorithm/pyqpanda_alg/QShadow/planner.py b/pyqpanda-algorithm/pyqpanda_alg/QShadow/planner.py new file mode 100644 index 0000000..ca245b5 --- /dev/null +++ b/pyqpanda-algorithm/pyqpanda_alg/QShadow/planner.py @@ -0,0 +1,321 @@ +"""Observable-aware local-Pauli measurement planning. + +The optimiser minimises a documented diagonal second-moment proxy, + + sum_j w_j c_j^2 / prod_{q in supp(P_j)} beta[q, P_j[q]], + +under independent single-qubit basis probabilities ``beta`` and a probability +floor. This objective follows the locally-biased classical-shadow motivation; +the block-coordinate solver in this module is an independent implementation. +It is a proxy, not a claim of globally optimal state-dependent variance. +""" + +from __future__ import annotations + +from collections import Counter +from dataclasses import dataclass +from typing import Sequence + +import numpy as np + +from .observable import PauliObservable, coerce_observables + +AXES = "XYZ" +_AXIS_TO_INDEX = {axis: index for index, axis in enumerate(AXES)} + + +def _validate_probabilities( + probabilities: Sequence[Sequence[float]] | np.ndarray, + *, + n_qubits: int | None = None, + probability_floor: float | None = None, +) -> np.ndarray: + result: np.ndarray = np.asarray(probabilities, dtype=float) + if result.ndim != 2 or result.shape[1] != 3: + raise ValueError("basis probabilities must have shape (n_qubits, 3)") + if n_qubits is not None and result.shape[0] != n_qubits: + raise ValueError( + f"probabilities describe {result.shape[0]} qubits, expected {n_qubits}" + ) + if not np.all(np.isfinite(result)): + raise ValueError("basis probabilities must be finite") + if np.any(result <= 0): + raise ValueError("every X/Y/Z basis probability must be positive") + if not np.allclose(result.sum(axis=1), 1.0, atol=1e-10): + raise ValueError("each qubit's X/Y/Z probabilities must sum to one") + if probability_floor is not None and np.any(result < probability_floor - 1e-12): + raise ValueError("initial probabilities violate probability_floor") + return result.copy() + + +def uniform_probabilities(n_qubits: int) -> np.ndarray: + """Return uniform local X/Y/Z probabilities for ``n_qubits``.""" + + if not isinstance(n_qubits, int) or isinstance(n_qubits, bool): + raise TypeError("n_qubits must be an integer") + if n_qubits <= 0: + raise ValueError("n_qubits must be positive") + return np.full((n_qubits, 3), 1.0 / 3.0, dtype=float) + + +def _weighted_terms( + observables: PauliObservable | Sequence[PauliObservable], + observable_weights: Sequence[float] | None, +) -> tuple[np.ndarray, np.ndarray, tuple[PauliObservable, ...]]: + checked = coerce_observables(observables) + weights: np.ndarray + if observable_weights is None: + weights = np.ones(len(checked), dtype=float) + else: + weights = np.asarray(observable_weights, dtype=float) + if weights.shape != (len(checked),): + raise ValueError("observable_weights must match the observable count") + if not np.all(np.isfinite(weights)) or np.any(weights < 0): + raise ValueError("observable_weights must be finite and non-negative") + if not np.any(weights > 0): + raise ValueError("at least one observable weight must be positive") + + axes: list[list[int]] = [] + term_weights: list[float] = [] + for observable, observable_weight in zip(checked, weights): + for term in observable.terms: + weight = float(observable_weight * term.coefficient**2) + if weight == 0.0 or term.is_identity: + continue + axes.append([_AXIS_TO_INDEX.get(axis, -1) for axis in term.pauli]) + term_weights.append(weight) + if not axes: + return ( + np.empty((0, checked[0].n_qubits), dtype=np.int8), + np.empty(0, dtype=float), + checked, + ) + return np.asarray(axes, dtype=np.int8), np.asarray(term_weights), checked + + +def _loss_from_arrays( + probabilities: np.ndarray, term_axes: np.ndarray, term_weights: np.ndarray +) -> float: + total = 0.0 + for axes, weight in zip(term_axes, term_weights): + coverage = 1.0 + for qubit, axis in enumerate(axes): + if axis >= 0: + coverage *= probabilities[qubit, axis] + total += weight / coverage + return float(total) + + +def planning_loss( + probabilities: Sequence[Sequence[float]], + observables: PauliObservable | Sequence[PauliObservable], + *, + observable_weights: Sequence[float] | None = None, +) -> float: + """Evaluate the diagonal second-moment proxy used by the optimiser.""" + + term_axes, term_weights, checked = _weighted_terms( + observables, observable_weights + ) + checked_probabilities = _validate_probabilities( + probabilities, n_qubits=checked[0].n_qubits + ) + return _loss_from_arrays(checked_probabilities, term_axes, term_weights) + + +@dataclass(frozen=True, slots=True) +class BasisOptimizationResult: + """Result and convergence diagnostics from basis-probability optimisation.""" + + probabilities: np.ndarray + initial_loss: float + final_loss: float + iterations: int + converged: bool + loss_history: tuple[float, ...] + probability_floor: float + + + +def _floor_constrained_inverse_cost(costs: np.ndarray, floor: float) -> np.ndarray: + """Solve ``min sum_a costs[a]/p[a]`` on a floor-constrained simplex.""" + + result: np.ndarray = np.full(3, floor, dtype=float) + free = [index for index, cost in enumerate(costs) if cost > 0.0] + if not free: + return np.full(3, 1.0 / 3.0) + + while free: + fixed_count = 3 - len(free) + remaining = 1.0 - fixed_count * floor + roots: np.ndarray = np.sqrt(costs[free]) + proposal = remaining * roots / roots.sum() + below_floor = [ + index for index, value in zip(free, proposal) if value < floor + ] + if not below_floor: + result[free] = proposal + break + free = [index for index in free if index not in below_floor] + + # Absorb round-off on the largest entry while preserving the floor. + result[int(np.argmax(result))] += 1.0 - float(result.sum()) + return result + + +def optimize_basis_probabilities( + observables: PauliObservable | Sequence[PauliObservable], + *, + observable_weights: Sequence[float] | None = None, + probability_floor: float = 0.02, + max_iterations: int = 200, + tolerance: float = 1e-10, + initial_probabilities: Sequence[Sequence[float]] | None = None, +) -> BasisOptimizationResult: + """Optimise independent local Pauli-basis probabilities. + + Each coordinate update exactly minimises the proxy with the other qubits + fixed. ``probability_floor`` keeps every basis observable and bounds inverse + propensities during post-processing. + """ + + if not 0.0 < probability_floor < 1.0 / 3.0: + raise ValueError("probability_floor must lie strictly between 0 and 1/3") + if not isinstance(max_iterations, int) or isinstance(max_iterations, bool): + raise TypeError("max_iterations must be an integer") + if max_iterations <= 0: + raise ValueError("max_iterations must be positive") + if not np.isfinite(tolerance) or tolerance <= 0: + raise ValueError("tolerance must be a finite positive number") + + term_axes, term_weights, checked = _weighted_terms( + observables, observable_weights + ) + n_qubits = checked[0].n_qubits + probabilities = ( + uniform_probabilities(n_qubits) + if initial_probabilities is None + else _validate_probabilities( + initial_probabilities, + n_qubits=n_qubits, + probability_floor=probability_floor, + ) + ) + initial_loss = _loss_from_arrays(probabilities, term_axes, term_weights) + history = [initial_loss] + if term_axes.shape[0] == 0: + return BasisOptimizationResult( + probabilities=probabilities, + initial_loss=0.0, + final_loss=0.0, + iterations=0, + converged=True, + loss_history=(0.0,), + probability_floor=probability_floor, + ) + + converged = False + iterations = 0 + for iteration in range(1, max_iterations + 1): + previous_loss = history[-1] + for qubit in range(n_qubits): + costs: np.ndarray = np.zeros(3, dtype=float) + for axes, weight in zip(term_axes, term_weights): + selected_axis = int(axes[qubit]) + if selected_axis < 0: + continue + other_coverage = 1.0 + for other_qubit, axis in enumerate(axes): + if other_qubit != qubit and axis >= 0: + other_coverage *= probabilities[other_qubit, axis] + costs[selected_axis] += weight / other_coverage + probabilities[qubit] = _floor_constrained_inverse_cost( + costs, probability_floor + ) + + current_loss = _loss_from_arrays(probabilities, term_axes, term_weights) + history.append(current_loss) + iterations = iteration + improvement = previous_loss - current_loss + if improvement < -1e-8 * max(1.0, abs(previous_loss)): + raise ArithmeticError("coordinate update unexpectedly increased planning loss") + if improvement <= tolerance * max(1.0, abs(previous_loss)): + converged = True + break + + return BasisOptimizationResult( + probabilities=probabilities.copy(), + initial_loss=initial_loss, + final_loss=history[-1], + iterations=iterations, + converged=converged, + loss_history=tuple(history), + probability_floor=probability_floor, + ) + + +@dataclass(frozen=True, slots=True) +class MeasurementPlan: + """A reproducible list of q0-first local-Pauli measurement bases. + + ``probabilities`` is one static ``(n_qubits, 3)`` distribution shared by + every shot in this plan. To use a different distribution in a later + round, create and execute another plan, then concatenate the resulting + :class:`~pyqpanda_alg.QShadow.dataset.ShadowDataset` objects. The datasets + retain the correct distribution for every shot. + """ + + bases: tuple[str, ...] + probabilities: np.ndarray + seed: int | None = None + + def __post_init__(self) -> None: + bases = tuple(self.bases) + if not bases: + raise ValueError("a measurement plan must contain at least one shot") + n_qubits = len(bases[0]) + if n_qubits == 0: + raise ValueError("basis strings cannot be empty") + for basis in bases: + if len(basis) != n_qubits or set(basis) - set(AXES): + raise ValueError("every basis must be an equal-length X/Y/Z string") + probabilities = _validate_probabilities( + self.probabilities, n_qubits=n_qubits + ) + object.__setattr__(self, "bases", bases) + object.__setattr__(self, "probabilities", probabilities) + + @property + def shots(self) -> int: + return len(self.bases) + + @property + def n_qubits(self) -> int: + return len(self.bases[0]) + + @property + def grouped_shots(self) -> dict[str, int]: + """Return identical basis strings compressed into backend jobs.""" + + return dict(Counter(self.bases)) + + +def sample_measurement_plan( + probabilities: Sequence[Sequence[float]], + shots: int, + *, + seed: int | None = None, +) -> MeasurementPlan: + """Sample a local-Pauli plan from independent per-qubit distributions.""" + + checked = _validate_probabilities(probabilities) + if not isinstance(shots, int) or isinstance(shots, bool): + raise TypeError("shots must be an integer") + if shots <= 0: + raise ValueError("shots must be positive") + rng = np.random.default_rng(seed) + sampled: np.ndarray = np.empty((shots, checked.shape[0]), dtype=np.int8) + for qubit in range(checked.shape[0]): + sampled[:, qubit] = rng.choice(3, size=shots, p=checked[qubit]) + bases = tuple("".join(AXES[index] for index in row) for row in sampled) + return MeasurementPlan(bases=bases, probabilities=checked, seed=seed) diff --git a/pyqpanda-algorithm/pyqpanda_alg/QShadow/sampler.py b/pyqpanda-algorithm/pyqpanda_alg/QShadow/sampler.py new file mode 100644 index 0000000..e6a7fce --- /dev/null +++ b/pyqpanda-algorithm/pyqpanda_alg/QShadow/sampler.py @@ -0,0 +1,261 @@ +"""Reference simulation and PyQPanda-compatible measurement execution.""" + +from __future__ import annotations + +from collections.abc import Callable, Mapping +from math import pi, sqrt +from typing import Any, Sequence + +import numpy as np + +from .dataset import ShadowDataset +from .planner import AXES, MeasurementPlan + +Counts = Mapping[str, int] +BasisRunner = Callable[[str, int], Counts] +ProgramExecutor = Callable[[Any, int], Counts] + +_HADAMARD = np.asarray([[1.0, 1.0], [1.0, -1.0]], dtype=complex) / sqrt(2.0) +_RX_PI_OVER_2 = np.asarray( + [[1.0, -1.0j], [-1.0j, 1.0]], dtype=complex +) / sqrt(2.0) +_IDENTITY: np.ndarray = np.eye(2, dtype=complex) + + +def _apply_single_qubit_gate( + state: np.ndarray, gate: np.ndarray, qubit: int +) -> np.ndarray: + result = state.copy() + stride = 1 << qubit + block = stride << 1 + for start in range(0, state.size, block): + for offset in range(stride): + zero = start + offset + one = zero + stride + amplitude_zero = state[zero] + amplitude_one = state[one] + result[zero] = gate[0, 0] * amplitude_zero + gate[0, 1] * amplitude_one + result[one] = gate[1, 0] * amplitude_zero + gate[1, 1] * amplitude_one + return result + + +def _basis_probabilities(state: np.ndarray, basis: str) -> np.ndarray: + rotated = state + for qubit, axis in enumerate(basis): + gate = _HADAMARD if axis == "X" else _RX_PI_OVER_2 if axis == "Y" else _IDENTITY + if axis != "Z": + rotated = _apply_single_qubit_gate(rotated, gate, qubit) + probabilities = np.abs(rotated) ** 2 + probabilities /= probabilities.sum() + return probabilities.real + + +def simulate_pauli_measurements( + statevector: Sequence[complex], + plan: MeasurementPlan, + *, + seed: int | None = None, +) -> ShadowDataset: + """Sample ``plan`` from an exact statevector using only NumPy. + + This independent reference path is useful for testing the estimator before + involving a simulator or QPU backend. It is exponential in qubit count and + intentionally intended only for small validation cases. + """ + + state: np.ndarray = np.asarray(statevector, dtype=np.complex128) + if state.ndim != 1 or state.size != 1 << plan.n_qubits: + raise ValueError( + f"statevector must have length {1 << plan.n_qubits} for this plan" + ) + norm = float(np.vdot(state, state).real) + if not np.all(np.isfinite(state)) or not np.isclose(norm, 1.0, atol=1e-10): + raise ValueError("statevector must be finite and normalised") + + rng = np.random.default_rng(seed) + probability_cache: dict[str, np.ndarray] = {} + bitstrings: list[str] = [] + for basis in plan.bases: + probabilities = probability_cache.get(basis) + if probabilities is None: + probabilities = _basis_probabilities(state, basis) + probability_cache[basis] = probabilities + outcome = int(rng.choice(state.size, p=probabilities)) + bitstrings.append( + "".join(str((outcome >> qubit) & 1) for qubit in range(plan.n_qubits)) + ) + return ShadowDataset.from_strings(plan.bases, bitstrings, plan.probabilities) + + +def _normalise_counts( + counts: Counts, + *, + n_qubits: int, + expected_shots: int, + bit_order: str, +) -> dict[str, int]: + if not isinstance(counts, Mapping): + raise TypeError("runner must return a mapping of bitstrings to integer counts") + if bit_order not in {"q0_first", "msb_first"}: + raise ValueError("bit_order must be 'q0_first' or 'msb_first'") + normalised: dict[str, int] = {} + for raw_bits, raw_count in counts.items(): + bits = str(raw_bits).replace(" ", "") + if len(bits) != n_qubits or set(bits) - {"0", "1"}: + raise ValueError(f"invalid {n_qubits}-qubit result key: {raw_bits!r}") + if not isinstance(raw_count, (int, np.integer)) or isinstance(raw_count, bool): + raise TypeError("measurement counts must be integers, not probabilities") + count = int(raw_count) + if count < 0: + raise ValueError("measurement counts cannot be negative") + if bit_order == "msb_first": + bits = bits[::-1] + normalised[bits] = normalised.get(bits, 0) + count + if sum(normalised.values()) != expected_shots: + raise ValueError( + f"runner returned {sum(normalised.values())} counts; " + f"expected {expected_shots}" + ) + return normalised + + +def run_measurement_plan( + plan: MeasurementPlan, + runner: BasisRunner, + *, + bit_order: str = "q0_first", + max_jobs: int | None = None, +) -> ShadowDataset: + """Execute each unique basis and convert counts into a QShadow dataset. + + Basis grouping may reorder shots within this call, which is safe because a + :class:`MeasurementPlan` has one static probability matrix shared by all of + its shots. Execute plans with different probability matrices separately + and combine their datasets with :meth:`ShadowDataset.concatenate`; this + preserves the appropriate propensity record for each round. + + ``max_jobs`` is a hard safety guard for remote or paid backends. No job is + submitted if the compressed plan exceeds the limit. + """ + + if not callable(runner): + raise TypeError("runner must be callable") + grouped = plan.grouped_shots + if max_jobs is not None: + if not isinstance(max_jobs, int) or isinstance(max_jobs, bool) or max_jobs <= 0: + raise ValueError("max_jobs must be a positive integer") + if len(grouped) > max_jobs: + raise RuntimeError( + f"plan needs {len(grouped)} backend jobs, exceeding max_jobs={max_jobs}" + ) + + bases: list[str] = [] + bitstrings: list[str] = [] + for basis, shots in grouped.items(): + counts = _normalise_counts( + runner(basis, shots), + n_qubits=plan.n_qubits, + expected_shots=shots, + bit_order=bit_order, + ) + for bits, count in counts.items(): + bases.extend([basis] * count) + bitstrings.extend([bits] * count) + return ShadowDataset.from_strings(bases, bitstrings, plan.probabilities) + + +def build_pyqpanda_program( + preparation: Callable[[Sequence[int]], Any], n_qubits: int, basis: str +) -> Any: + """Build a measured PyQPanda3 program for one q0-first basis string.""" + + if not callable(preparation): + raise TypeError("preparation must be callable") + if not isinstance(n_qubits, int) or isinstance(n_qubits, bool) or n_qubits <= 0: + raise ValueError("n_qubits must be a positive integer") + if len(basis) != n_qubits or set(basis) - set(AXES): + raise ValueError("basis must be an n_qubits-long X/Y/Z string") + + # Lazy import keeps pure statistical post-processing independently testable. + from pyqpanda3.core import H, RX, QProg, measure + + program = QProg(n_qubits) + qubits = program.qubits() + prepared = preparation(qubits) + if prepared is not None: + program << prepared + for qubit, axis in zip(qubits, basis): + if axis == "X": + program << H(qubit) + elif axis == "Y": + # RX(+pi/2)^dagger Z RX(+pi/2) = Y for pre-measurement rotation. + program << RX(qubit, pi / 2.0) + for qubit in qubits: + program << measure(qubit, qubit) + return program + + +class PyQPandaRunner: + """Create basis programs and execute them locally or through a callback. + + Parameters + ---------- + preparation: + Callable receiving the program's qubit indices and returning a QCircuit, + QProg, gate, or ``None``. Returning ``None`` prepares ``|0...0>``. + n_qubits: + Program width. + executor: + Optional ``executor(program, shots) -> counts`` callback. If omitted, + a local :class:`pyqpanda3.core.CPUQVM` is used. A cloud adapter should + capture an already-authenticated backend here; QShadow never accepts or + stores API keys. + backend_bit_order: + PyQPanda count dictionaries are normally MSB-first. Results are + converted to QShadow's q0-first convention before return. + """ + + def __init__( + self, + preparation: Callable[[Sequence[int]], Any], + n_qubits: int, + *, + executor: ProgramExecutor | None = None, + backend_bit_order: str = "msb_first", + ) -> None: + if not callable(preparation): + raise TypeError("preparation must be callable") + if not isinstance(n_qubits, int) or isinstance(n_qubits, bool) or n_qubits <= 0: + raise ValueError("n_qubits must be a positive integer") + if backend_bit_order not in {"q0_first", "msb_first"}: + raise ValueError("backend_bit_order must be 'q0_first' or 'msb_first'") + if executor is not None and not callable(executor): + raise TypeError("executor must be callable") + self.preparation = preparation + self.n_qubits = n_qubits + self.backend_bit_order = backend_bit_order + self._machine: Any | None = None + if executor is None: + from pyqpanda3.core import CPUQVM + + self._machine = CPUQVM() + self.executor: ProgramExecutor = self._execute_local + else: + self.executor = executor + + def _execute_local(self, program: Any, shots: int) -> Counts: + assert self._machine is not None + self._machine.run(program, shots=shots) + return self._machine.result().get_counts() + + def __call__(self, basis: str, shots: int) -> dict[str, int]: + if not isinstance(shots, int) or isinstance(shots, bool) or shots <= 0: + raise ValueError("shots must be a positive integer") + program = build_pyqpanda_program(self.preparation, self.n_qubits, basis) + counts = self.executor(program, shots) + return _normalise_counts( + counts, + n_qubits=self.n_qubits, + expected_shots=shots, + bit_order=self.backend_bit_order, + ) diff --git a/pyqpanda-algorithm/pyqpanda_alg/__init__.py b/pyqpanda-algorithm/pyqpanda_alg/__init__.py index 12d6808..33a33b3 100644 --- a/pyqpanda-algorithm/pyqpanda_alg/__init__.py +++ b/pyqpanda-algorithm/pyqpanda_alg/__init__.py @@ -46,6 +46,7 @@ from . import QUBO from . import QCmp from . import QAE +from . import QShadow from . import QSVD from . import QSVR from . import Grover diff --git a/test/QShadow/Test_estimator.py b/test/QShadow/Test_estimator.py new file mode 100644 index 0000000..4739e98 --- /dev/null +++ b/test/QShadow/Test_estimator.py @@ -0,0 +1,133 @@ +import numpy as np +import pytest + +from pyqpanda_alg.QShadow import ( + PauliObservable, + PauliTerm, + ShadowBudgetController, + ShadowDataset, + empirical_bernstein_radius, + estimate_observable, + estimate_observables, + optimize_basis_probabilities, + sample_measurement_plan, + simulate_pauli_measurements, + term_shadow_samples, + uniform_probabilities, +) + + +def test_term_samples_use_matches_signs_and_inverse_propensities(): + probabilities = np.asarray([[0.2, 0.3, 0.5]]) + dataset = ShadowDataset.from_strings( + ["X", "Y", "Z", "Z"], ["0", "0", "0", "1"], probabilities + ) + samples = term_shadow_samples(dataset, PauliTerm("Z", 2.0)) + assert np.allclose(samples, [0.0, 0.0, 4.0, -4.0]) + + +def test_identity_observable_is_exact_for_every_dataset(): + plan = sample_measurement_plan(uniform_probabilities(2), 30, seed=1) + dataset = simulate_pauli_measurements([1, 0, 0, 0], plan, seed=2) + identity = PauliObservable.from_terms({"II": 2.5}, name="identity") + estimate = estimate_observable(dataset, identity) + assert estimate.value == 2.5 + assert estimate.half_width == 0.0 + assert estimate.minimum_term_matches == 30 + + +def test_monte_carlo_estimate_is_unbiased_with_fixed_seed_regression(): + observable = PauliObservable.from_terms( + {"ZI": 0.8, "IZ": -0.4, "XX": 0.2, "II": 0.1}, name="hamiltonian" + ) + state = np.asarray([1, 0, 0, 0], dtype=complex) + probabilities = optimize_basis_probabilities(observable).probabilities + plan = sample_measurement_plan(probabilities, 12_000, seed=10) + dataset = simulate_pauli_measurements(state, plan, seed=11) + estimate = estimate_observable(dataset, observable, confidence=0.95) + exact = observable.exact_expectation(state) + assert estimate.value == pytest.approx(exact, abs=0.035) + assert estimate.lower <= exact <= estimate.upper + + +def test_concatenated_rounds_preserve_each_round_probabilities(): + observable = PauliObservable.from_terms({"Z": 1.0}, name="z") + first_plan = sample_measurement_plan(uniform_probabilities(1), 4_000, seed=20) + second_probabilities = np.asarray([[0.05, 0.05, 0.90]]) + second_plan = sample_measurement_plan(second_probabilities, 4_000, seed=21) + state = [1, 0] + combined = ShadowDataset.concatenate( + simulate_pauli_measurements(state, first_plan, seed=22), + simulate_pauli_measurements(state, second_plan, seed=23), + ) + estimate = estimate_observable(combined, observable) + assert combined.probabilities[0, 0, 2] == pytest.approx(1 / 3) + assert combined.probabilities[-1, 0, 2] == pytest.approx(0.90) + assert estimate.value == pytest.approx(1.0, abs=0.04) + + +def test_many_observables_use_simultaneous_bonferroni_delta(): + plan = sample_measurement_plan(uniform_probabilities(1), 200, seed=30) + dataset = simulate_pauli_measurements([1, 0], plan, seed=31) + observables = [ + PauliObservable.from_terms({"Z": 1}, name="z"), + PauliObservable.from_terms({"X": 1}, name="x"), + ] + estimates = estimate_observables(dataset, observables, family_confidence=0.90) + assert set(estimates) == {"z", "x"} + assert estimates["z"].failure_probability == pytest.approx(0.05) + assert estimates["x"].failure_probability == pytest.approx(0.05) + + +def test_empirical_bernstein_radius_shrinks_with_repeated_data(): + short = empirical_bernstein_radius( + [1.0, -1.0] * 10, value_range=2.0, failure_probability=0.05 + ) + long = empirical_bernstein_radius( + [1.0, -1.0] * 1000, value_range=2.0, failure_probability=0.05 + ) + assert long < short + + +def test_budget_controller_completes_exact_identity_target(): + plan = sample_measurement_plan(uniform_probabilities(1), 100, seed=40) + dataset = simulate_pauli_measurements([1, 0], plan, seed=41) + identity = PauliObservable.from_terms({"I": 1}, name="identity") + report = ShadowBudgetController( + identity, target_half_width=1e-6, min_shots=100, max_shots=500 + ).evaluate(dataset) + assert report.complete + assert not report.exhausted + assert report.next_batch_shots == 0 + assert report.reason == "target_reached" + + +def test_budget_controller_never_exceeds_hard_maximum(): + plan = sample_measurement_plan(uniform_probabilities(1), 100, seed=50) + dataset = simulate_pauli_measurements([1, 0], plan, seed=51) + observable = PauliObservable.from_terms({"Z": 1}, name="z") + report = ShadowBudgetController( + observable, + target_half_width=1e-9, + min_shots=50, + max_shots=150, + batch_shots=200, + ).evaluate(dataset) + assert not report.complete + assert report.next_batch_shots == 50 + + exhausted_data = ShadowDataset.concatenate( + dataset, + simulate_pauli_measurements( + [1, 0], sample_measurement_plan(uniform_probabilities(1), 50, seed=52), seed=53 + ), + ) + exhausted = ShadowBudgetController( + observable, + target_half_width=1e-9, + min_shots=50, + max_shots=150, + ).evaluate(exhausted_data) + assert exhausted.exhausted + assert exhausted.next_batch_shots == 0 + assert exhausted.reason == "max_shots_reached" diff --git a/test/QShadow/Test_observable.py b/test/QShadow/Test_observable.py new file mode 100644 index 0000000..860c185 --- /dev/null +++ b/test/QShadow/Test_observable.py @@ -0,0 +1,58 @@ +import numpy as np +import pytest + +from pyqpanda_alg.QShadow import PauliObservable, PauliTerm, pauli_expectation + + +def test_pauli_term_normalises_and_reports_support(): + term = PauliTerm(" x i z ", 2) + assert term.pauli == "XIZ" + assert term.support == (0, 2) + assert term.weight == 2 + assert term.n_qubits == 3 + + +def test_observable_combines_duplicate_terms(): + observable = PauliObservable.from_terms( + [("ZI", 0.5), PauliTerm("ZI", 1.25), ("IZ", -0.25)], name="energy" + ) + assert observable.name == "energy" + assert observable.terms == (PauliTerm("ZI", 1.75), PauliTerm("IZ", -0.25)) + + +def test_pauli_term_rejects_invalid_or_non_hermitian_input(): + with pytest.raises(ValueError, match="invalid Pauli"): + PauliTerm("XA") + with pytest.raises(ValueError, match="must be real"): + PauliTerm("X", 1 + 2j) + with pytest.raises(ValueError, match="finite"): + PauliTerm("X", np.inf) + + +def test_observable_rejects_mixed_width_and_duplicate_names_are_handled_elsewhere(): + with pytest.raises(ValueError, match="equal length"): + PauliObservable((PauliTerm("X"), PauliTerm("ZZ"))) + + +def test_exact_single_qubit_pauli_expectations(): + plus = np.asarray([1.0, 1.0]) / np.sqrt(2.0) + plus_i = np.asarray([1.0, 1.0j]) / np.sqrt(2.0) + assert pauli_expectation([1, 0], "Z") == pytest.approx(1.0) + assert pauli_expectation(plus, "X") == pytest.approx(1.0) + assert pauli_expectation(plus_i, "Y") == pytest.approx(1.0) + + +def test_exact_expectation_uses_q0_first_pauli_strings(): + # Statevector index 1 is |q1 q0> = |01>: q0=1 and q1=0. + state = np.asarray([0, 1, 0, 0], dtype=complex) + observable = PauliObservable.from_terms( + {"ZI": 2.0, "IZ": 0.5, "II": 0.25}, name="q0-first" + ) + assert observable.exact_expectation(state) == pytest.approx(-1.25) + + +def test_exact_expectation_rejects_wrong_or_unnormalised_statevector(): + with pytest.raises(ValueError, match="expected 4"): + pauli_expectation([1, 0], "ZZ") + with pytest.raises(ValueError, match="normalised"): + pauli_expectation([2, 0], "Z") diff --git a/test/QShadow/Test_planner.py b/test/QShadow/Test_planner.py new file mode 100644 index 0000000..d6a3e31 --- /dev/null +++ b/test/QShadow/Test_planner.py @@ -0,0 +1,95 @@ +import numpy as np +import pytest + +from pyqpanda_alg.QShadow import ( + MeasurementPlan, + PauliObservable, + optimize_basis_probabilities, + planning_loss, + sample_measurement_plan, + uniform_probabilities, +) + + +def test_uniform_probabilities_shape_and_values(): + probabilities = uniform_probabilities(4) + assert probabilities.shape == (4, 3) + assert np.allclose(probabilities, 1 / 3) + + +def test_single_axis_optimum_saturates_unused_axes_at_floor(): + observable = PauliObservable.from_terms({"Z": 1.0}, name="z") + result = optimize_basis_probabilities(observable, probability_floor=0.05) + assert result.converged + assert np.allclose(result.probabilities[0], [0.05, 0.05, 0.90]) + assert result.final_loss < result.initial_loss + + +def test_unused_qubit_stays_uniform(): + observable = PauliObservable.from_terms({"ZI": 1.0}, name="local-z") + result = optimize_basis_probabilities(observable, probability_floor=0.02) + assert np.allclose(result.probabilities[1], [1 / 3, 1 / 3, 1 / 3]) + + +def test_optimizer_reduces_multi_observable_proxy_and_respects_floor(): + observables = [ + PauliObservable.from_terms( + {"ZZI": 1.4, "IZZ": -0.7, "XII": 0.2}, name="energy" + ), + PauliObservable.from_terms( + {"ZIZ": 1.0, "YYY": 0.1}, name="correlator" + ), + ] + uniform = uniform_probabilities(3) + result = optimize_basis_probabilities( + observables, + observable_weights=[2.0, 0.5], + probability_floor=0.03, + ) + assert np.all(result.probabilities >= 0.03 - 1e-12) + assert np.allclose(result.probabilities.sum(axis=1), 1.0) + assert result.final_loss == pytest.approx( + planning_loss( + result.probabilities, observables, observable_weights=[2.0, 0.5] + ) + ) + assert result.final_loss < planning_loss( + uniform, observables, observable_weights=[2.0, 0.5] + ) + assert all( + later <= earlier + 1e-10 + for earlier, later in zip(result.loss_history, result.loss_history[1:]) + ) + + +def test_zero_observable_returns_uniform_without_iteration(): + observable = PauliObservable.from_terms({"II": 0.0}, name="zero") + result = optimize_basis_probabilities(observable) + assert result.iterations == 0 + assert result.final_loss == 0.0 + assert np.allclose(result.probabilities, 1 / 3) + + +def test_sample_plan_is_reproducible_and_grouped_counts_sum_to_shots(): + probabilities = np.asarray([[0.2, 0.3, 0.5], [0.6, 0.2, 0.2]]) + first = sample_measurement_plan(probabilities, 200, seed=123) + second = sample_measurement_plan(probabilities, 200, seed=123) + assert first.bases == second.bases + assert first.shots == 200 + assert first.n_qubits == 2 + assert sum(first.grouped_shots.values()) == 200 + + +def test_measurement_plan_rejects_inconsistent_input(): + with pytest.raises(ValueError, match="equal-length"): + MeasurementPlan(("XZ", "X"), uniform_probabilities(2)) + with pytest.raises(ValueError, match="sum to one"): + MeasurementPlan(("X",), np.asarray([[0.2, 0.2, 0.2]])) + + +def test_optimizer_rejects_invalid_floor_and_weights(): + observable = PauliObservable.from_terms({"X": 1}, name="x") + with pytest.raises(ValueError, match="between 0 and 1/3"): + optimize_basis_probabilities(observable, probability_floor=1 / 3) + with pytest.raises(ValueError, match="match"): + optimize_basis_probabilities(observable, observable_weights=[1, 2]) diff --git a/test/QShadow/Test_sampler.py b/test/QShadow/Test_sampler.py new file mode 100644 index 0000000..d8c12be --- /dev/null +++ b/test/QShadow/Test_sampler.py @@ -0,0 +1,118 @@ +import numpy as np +import pytest + +from pyqpanda3.core import H, QCircuit, RX, X + +from pyqpanda_alg.QShadow import ( + MeasurementPlan, + PyQPandaRunner, + ShadowDataset, + build_pyqpanda_program, + run_measurement_plan, + simulate_pauli_measurements, + uniform_probabilities, +) + + +def test_reference_sampler_gets_deterministic_x_and_y_eigenvalues(): + plus = np.asarray([1.0, 1.0]) / np.sqrt(2.0) + plus_i = np.asarray([1.0, 1.0j]) / np.sqrt(2.0) + x_plan = MeasurementPlan(("X",) * 50, uniform_probabilities(1)) + y_plan = MeasurementPlan(("Y",) * 50, uniform_probabilities(1)) + assert set(simulate_pauli_measurements(plus, x_plan, seed=1).outcomes[:, 0]) == {1} + assert set(simulate_pauli_measurements(plus_i, y_plan, seed=2).outcomes[:, 0]) == {1} + + +def test_fake_runner_q0_first_and_msb_first_conversion(): + plan = MeasurementPlan(("ZZ",), uniform_probabilities(2)) + q0_first = run_measurement_plan(plan, lambda basis, shots: {"01": 1}) + msb_first = run_measurement_plan( + plan, lambda basis, shots: {"01": 1}, bit_order="msb_first" + ) + assert q0_first.bitstrings == ("01",) + assert msb_first.bitstrings == ("10",) + + +def test_runner_rejects_count_mismatch_and_non_integer_counts(): + plan = MeasurementPlan(("Z", "Z"), uniform_probabilities(1)) + with pytest.raises(ValueError, match="expected 2"): + run_measurement_plan(plan, lambda basis, shots: {"0": 1}) + with pytest.raises(TypeError, match="integers"): + run_measurement_plan(plan, lambda basis, shots: {"0": 1.0, "1": 1.0}) + + +def test_max_jobs_guard_runs_before_any_remote_callback(): + plan = MeasurementPlan(("X", "Y", "Z"), uniform_probabilities(1)) + called = False + + def runner(basis, shots): + nonlocal called + called = True + return {"0": shots} + + with pytest.raises(RuntimeError, match="exceeding max_jobs"): + run_measurement_plan(plan, runner, max_jobs=2) + assert not called + + +def test_grouped_remote_rounds_preserve_each_round_probabilities(): + first_probabilities = np.asarray([[0.2, 0.3, 0.5]]) + second_probabilities = np.asarray([[0.7, 0.2, 0.1]]) + first = run_measurement_plan( + MeasurementPlan(("Z", "X", "Z"), first_probabilities), + lambda basis, shots: {"0": shots}, + ) + second = run_measurement_plan( + MeasurementPlan(("X", "Z", "Z"), second_probabilities), + lambda basis, shots: {"0": shots}, + ) + combined = ShadowDataset.concatenate(first, second) + + selected_propensities = combined.probabilities[ + np.arange(combined.shots), 0, combined.bases[:, 0] + ] + assert combined.basis_strings == ("Z", "Z", "X", "X", "Z", "Z") + assert np.allclose(selected_propensities, [0.5, 0.5, 0.2, 0.7, 0.1, 0.1]) + + +def test_cpuqvm_runner_converts_pyqpanda_msb_order_to_q0_first(): + plan = MeasurementPlan(("ZZ",) * 20, uniform_probabilities(2)) + runner = PyQPandaRunner(lambda qubits: QCircuit() << X(qubits[0]), 2) + dataset = run_measurement_plan(plan, runner, max_jobs=1) + # q0=1 and q1=0. Native PyQPanda counts are '01'; runner returns q0-first '10'. + assert set(dataset.bitstrings) == {"10"} + + +def test_cpuqvm_runner_x_and_y_basis_signs(): + x_plan = MeasurementPlan(("X",) * 30, uniform_probabilities(1)) + y_plan = MeasurementPlan(("Y",) * 30, uniform_probabilities(1)) + x_runner = PyQPandaRunner(lambda qubits: QCircuit() << H(qubits[0]), 1) + y_runner = PyQPandaRunner( + lambda qubits: QCircuit() << RX(qubits[0], -np.pi / 2), 1 + ) + assert set(run_measurement_plan(x_plan, x_runner).outcomes[:, 0]) == {1} + assert set(run_measurement_plan(y_plan, y_runner).outcomes[:, 0]) == {1} + + +def test_custom_executor_receives_program_and_shots_without_credentials(): + calls = [] + + def executor(program, shots): + calls.append((type(program).__name__, shots)) + return {"0": shots} + + plan = MeasurementPlan(("Z",) * 7, uniform_probabilities(1)) + runner = PyQPandaRunner( + lambda qubits: None, + 1, + executor=executor, + backend_bit_order="msb_first", + ) + dataset = run_measurement_plan(plan, runner, max_jobs=1) + assert calls == [("QProg", 7)] + assert dataset.bitstrings == ("0",) * 7 + + +def test_pyqpanda_builder_validates_basis_width(): + with pytest.raises(ValueError, match="X/Y/Z"): + build_pyqpanda_program(lambda qubits: None, 2, "X") diff --git a/test/pytest.ini b/test/pytest.ini index 2887f64..d4f2f48 100644 --- a/test/pytest.ini +++ b/test/pytest.ini @@ -5,6 +5,7 @@ testpaths = QRAM QPCA QSVM + QShadow python_files =