From 076f7a1c1b4f28553bb9674376db655876ec792c Mon Sep 17 00:00:00 2001 From: Ruge Lin Date: Tue, 21 Jul 2026 15:00:13 +0800 Subject: [PATCH] Add exact-logical QBP validation --- PAPER1_NUMERICS.md | 225 ++++++++ README.md | 503 ++---------------- qbp_validation/README.md | 104 ++++ qbp_validation/__init__.py | 3 + qbp_validation/cases.py | 111 ++++ qbp_validation/circuits.py | 351 ++++++++++++ qbp_validation/conventions.py | 212 ++++++++ qbp_validation/decoders.py | 115 ++++ qbp_validation/reference.py | 280 ++++++++++ qbp_validation/tests/__init__.py | 1 + qbp_validation/tests/test_checkpoints.py | 59 ++ .../tests/test_complex_magnitude.py | 60 +++ qbp_validation/tests/test_complex_phase.py | 47 ++ qbp_validation/tests/test_conventions.py | 54 ++ qbp_validation/tests/test_decoders.py | 57 ++ .../tests/test_four_qubit_example.py | 60 +++ qbp_validation/tests/test_polyspherical.py | 77 +++ qbp_validation/tests/test_real_frame.py | 50 ++ .../tests/test_real_global_estimator.py | 43 ++ qbp_validation/tests/test_resource_ledger.py | 32 ++ qbp_validation/tests/test_singular_cases.py | 60 +++ 21 files changed, 2036 insertions(+), 468 deletions(-) create mode 100644 PAPER1_NUMERICS.md create mode 100644 qbp_validation/README.md create mode 100644 qbp_validation/__init__.py create mode 100644 qbp_validation/cases.py create mode 100644 qbp_validation/circuits.py create mode 100644 qbp_validation/conventions.py create mode 100644 qbp_validation/decoders.py create mode 100644 qbp_validation/reference.py create mode 100644 qbp_validation/tests/__init__.py create mode 100644 qbp_validation/tests/test_checkpoints.py create mode 100644 qbp_validation/tests/test_complex_magnitude.py create mode 100644 qbp_validation/tests/test_complex_phase.py create mode 100644 qbp_validation/tests/test_conventions.py create mode 100644 qbp_validation/tests/test_decoders.py create mode 100644 qbp_validation/tests/test_four_qubit_example.py create mode 100644 qbp_validation/tests/test_polyspherical.py create mode 100644 qbp_validation/tests/test_real_frame.py create mode 100644 qbp_validation/tests/test_real_global_estimator.py create mode 100644 qbp_validation/tests/test_resource_ledger.py create mode 100644 qbp_validation/tests/test_singular_cases.py diff --git a/PAPER1_NUMERICS.md b/PAPER1_NUMERICS.md new file mode 100644 index 0000000..a45ee5b --- /dev/null +++ b/PAPER1_NUMERICS.md @@ -0,0 +1,225 @@ +# Paper I numerical material + +This page documents the repository material for +**[A Compass on the Quantum State Sphere: The Hopf Ansatz for Arbitrary Pure-State Optimization](https://arxiv.org/abs/2607.14231)**. + +The first paper studies the Hopf ansatz as an optimization chart. The repository contains +real- and complex-Hopf stress tests, adaptive baselines, diagnostics, plotting scripts, +finite-shot estimator checks, CNOT-ledger checks, and a small Qibo VQE circuit demo. +The generated CSV archives are derived artifacts and are intentionally not tracked. + +## Main files + +| File | Role | +|---|---| +| `hopf_utils.py` | Real and complex Hopf maps, inverse maps, Jacobians, diagonal metrics, tangent-state assignments, gate schedules, and optional Qibo checks. | +| `hopf_data.py` | Generates the multi-size real-Hopf synthetic stress-test datasets. | +| `adam_data.py` | Generates the adaptive Adam baseline datasets. | +| `diagnose_hopf.py` | Checks completeness and numerical quality of the Hopf datasets. | +| `diagnose_adam.py` | Checks completeness and numerical quality of the Adam datasets. | +| `plot_hopf.py` | Generates the real-Hopf and Adam summary figures. | +| `hopf_complex.py` | Focused complex-Hopf stress test at `n=6`. | +| `finite_shot_sanity_check.py` | Fixed-state finite-shot check of the symmetric signed-branch gradient identity. | +| `hopf_gate_count.py` | Checks generated real and complex gate schedules against the closed-form CNOT ledger. | +| `VQE_qibo.py` | Local `n=4` real/complex Qibo and statevector layerwise-gradient demonstration. | + +Included small release figures are regenerated by their corresponding scripts: + +```text +hopf_complex.png +finite_shot_sanity.png +VQE_qibo.png +``` + +## Scope + +The synthetic tasks comprise three VQE families and three metrology-inspired families: + +- parent-Hamiltonian recovery; +- hidden Hamming-spectrum minimization; +- small-gap spectrum minimization; +- single-target Fisher proxy; +- normalized pure-state QFI; +- balanced two-target Fisher proxy. + +These experiments test optimization behavior. They are separate from the second +manuscript's deterministic theorem-validation suite in `qbp_validation/`. + +## Requirements + +Use Python 3.10 or newer. + +```bash +python -m pip install -r requirements.txt +``` + +Explicit Qibo circuit paths additionally use: + +```bash +python -m pip install -r requirements-optional.txt +``` + +## Minimal safeguards + +```bash +python hopf_utils.py + +mkdir -p figures_safeguards +python hopf_gate_count.py \ + --nmin 2 \ + --nmax 6 \ + --out figures_safeguards/hopf_gate_count_smoke.pdf + +MPLBACKEND=Agg python VQE_qibo.py \ + --steps 2 \ + --shots 20 \ + --sampler statevector \ + --output figures_safeguards/VQE_qibo_smoke.png \ + --log-every 0 + +python hopf_complex.py \ + --quick \ + --outdir figures_safeguards \ + --csv-name complex_hopf_smoke.csv \ + --plot-name hopf_complex_smoke.png + +python finite_shot_sanity_check.py \ + --shots 100 1000 \ + --trials 3 \ + --output-dir figures_safeguards +``` + +Use `--sampler qibo-explicit` in `VQE_qibo.py` to require the explicit Qibo path. +The default `auto` mode uses Qibo when available and otherwise falls back to the +statevector implementation. + +## Full real-Hopf datasets + +```bash +mkdir -p data + +for n in 6 7 8 9 10; do + python hopf_data.py --n "$n" --steps 200 --num-seeds 10 --outdir data + python adam_data.py --n "$n" --steps 200 --num-seeds 10 --outdir data +done +``` + +The geometry-native tracks are: + +```text +Hopf-EGT-CG +Hopf-Riemannian-LBFGS +Hopf-Riemannian-BB +``` + +The adaptive baselines are: + +```text +Hopf-Adam +Mottonen-ideal-PS-Adam +``` + +## Diagnostics + +```bash +repo_root="$(pwd)" +mkdir -p diagnostics + +python diagnose_hopf.py \ + --indir data \ + --ns 6-10 \ + --steps 200 \ + --num-seeds 10 \ + --out "$repo_root/diagnostics/hopf_data_diagnostics.txt" + +python diagnose_adam.py \ + --indir data \ + --ns 6-10 \ + --steps 200 \ + --num-seeds 10 \ + --out "$repo_root/diagnostics/adam_data_diagnostics.txt" +``` + +The reports check expected files and rows, task/mode/seed/step coverage, parse errors, +vector lengths, nonfinite values, state normalization, convergence summaries, and final +rankings. The diagnostic scripts stream the large vector columns rather than loading the +full archive into pandas. + +## Figures + +```bash +mkdir -p figures +python plot_hopf.py \ + --indir data \ + --adam-dir data \ + --outdir figures \ + --ns 6-10 \ + --detail-n 10 \ + --formats pdf png +``` + +## Focused complex-Hopf run + +```bash +python hopf_complex.py +``` + +This generates: + +```text +complex_hopf_stress_data.csv +hopf_complex.png +``` + +For a smaller run: + +```bash +python hopf_complex.py --quick +``` + +For plotting only from an existing CSV: + +```bash +python hopf_complex.py --plot-only +``` + +## Standalone safeguards + +### CNOT ledger + +```bash +python hopf_gate_count.py +``` + +This compares generated schedules with the closed-form real and complex CNOT-count +formulas. The resulting plot is a local safeguard, not a hardware compilation benchmark. + +### Finite-shot signed-branch identity + +```bash +python finite_shot_sanity_check.py +``` + +The script fixes one `n=6` real Hopf state, verifies the exact branch identity against an +independent tree-gradient reference, and then samples the two branch energies at several +shot counts. This is an estimator sanity check, not a total full-gradient sampling study. + +### Local Qibo VQE demonstration + +```bash +MPLBACKEND=Agg python VQE_qibo.py +``` + +The script runs local `n=4` real and complex VQE toys and compares exact Hopf-gradient +Adam with sampled layerwise-circuit Adam. It is a local circuit-realizability demonstration, +not the global-frame implementation of the second manuscript and not an asymptotic gate-count +benchmark. + +## Determinism and complete commands + +The scripts use deterministic seeds by default. Exact regeneration commands, seed rules, +chart tolerances, and output conventions are listed in +[`REPRODUCIBILITY.md`](REPRODUCIBILITY.md). + +The repository landing page before the two-paper reorganization is preserved in Git history +at commit `1188e8a1af5326d62d8c788499c1d59a06fe450c`. diff --git a/README.md b/README.md index 34b1495..639b56d 100644 --- a/README.md +++ b/README.md @@ -1,491 +1,58 @@ -# Hopf Ansatz Stress-Test Code +# Hopf Ansatz -This repository contains reproducible synthetic stress-test code for the Hopf ansatz. +This repository accompanies two papers built on the same Hopf-coordinate construction. +The numerical optimization studies of the first paper and the exact-logical validation +of the second paper are kept in one repository but are separated by scope. -The code generates real- and complex-Hopf VQE and metrology-inspired optimization traces, runs seed-aware diagnostics, produces GitHub-facing comparison plots, and includes standalone safeguards for the Hopf CNOT-count formulas, the finite-shot signed-branch gradient estimator, and Qibo layerwise gradient-access circuits. +| Paper | Scope | Repository material | +|---|---|---| +| **[A Compass on the Quantum State Sphere: The Hopf Ansatz for Arbitrary Pure-State Optimization](https://arxiv.org/abs/2607.14231)** | Hopf coordinates, inverse maps, geometry, optimization algorithms, synthetic VQE and metrology-inspired stress tests | Existing root-level scripts, diagnostics, plots, and [`PAPER1_NUMERICS.md`](PAPER1_NUMERICS.md) | +| **The Compass in Reverse: Quantum Backpropagation with the Hopf Ansatz** | Exact-logical complete-gradient estimators, global frames, complex phase decoding, polyspherical frames, and checkpointed adjoints | Deterministic Qibo circuit and decoder checks in [`qbp_validation/`](qbp_validation/) | -Full generated CSV datasets are intentionally not tracked in GitHub because they are derived artifacts. They can be regenerated from deterministic scripts; see `REPRODUCIBILITY.md`. The included `VQE_qibo.png`, `hopf_complex.png`, and `finite_shot_sanity.png` files are small release artifacts that can be regenerated from `VQE_qibo.py`, `hopf_complex.py`, and `finite_shot_sanity_check.py`, respectively. Diagnostic text files and generated CSV files are ignored by default. +## Scope separation -## Repository contents +The first paper contains numerical optimization studies. Its scripts generate optimizer +traces, diagnostics, and figures. -| File | Role | -|---|---| -| `hopf_utils.py` | Core Hopf utilities: real and complex coordinate maps, inverse maps, Jacobians, diagonal metrics, tangent-state assignments, gate schedules, and optional Qibo circuit checks. | -| `hopf_data.py` | Generates synthetic stress-test CSVs for the three geometry-native Hopf optimizers. | -| `adam_data.py` | Generates synthetic stress-test CSVs for the two adaptive Adam baselines. | -| `diagnose_hopf.py` | Checks completeness and numerical quality of Hopf optimizer CSVs. | -| `diagnose_adam.py` | Checks completeness and numerical quality of Adam baseline CSVs. | -| `plot_hopf.py` | Generates summary and convergence plots from the real-Hopf and Adam CSVs. | -| `hopf_complex.py` | Focused complex-Hopf stress test at `n=6`. It runs the same six VQE and metrology-inspired objective families with complex scrambling, ten deterministic initial states per task, and four optimizer tracks; it writes the CSV, prints diagnostics, and produces the final-gap summary plot. | -| `hopf_complex.png` | Included final-gap summary for the default complex-Hopf stress test. | -| `finite_shot_sanity_check.py` | Fixed-state finite-shot check of the symmetric signed-branch Hopf gradient estimator. It compares sampled full coordinate gradients with an exact tree-gradient reference, prints statistical diagnostics, and generates its CSV and requested plot formats. | -| `finite_shot_sanity.png` | Included finite-shot estimator summary for the default `n=6`, 50-trial experiment. | -| `hopf_gate_count.py` | Safeguard script for the CNOT-count formulas. It builds real and complex Hopf gate schedules and checks numerical counts against the closed-form binomial formulas. | -| `VQE_qibo.py` | Functional Qibo/statevector layerwise gradient-access safeguard for local `n=4` real and complex Hopf VQE toys; compares exact Hopf-gradient Adam with sampled layerwise-circuit Adam. This is a local circuit-realizability demo, not the asymptotically optimized indexed-gradient scaling implementation. | -| `VQE_qibo.png` | Included output panel from `VQE_qibo.py`, showing the real and complex `n=4` local toy trajectories. | +The second paper is theoretical and makes no numerical-performance claim. Its repository +material consists only of deterministic sanity checks of the exact-logical state, frame, +measurement, decoder, singular-chart, checkpoint, and compiler-ledger identities. It does +not contain optimization comparisons, empirical scaling plots, hardware experiments, or +noise studies. -## Scope relative to the paper +## Repository map -The paper is kept focused on the analytic construction. This repository carries implementation-level material: inverse-map checks, metric/Jacobian consistency checks, tangent-state synthesis checks, generated gate schedules, CNOT-count safeguards, the multi-size real-Hopf optimizer traces, a focused complex-Hopf stress test, a finite-shot signed-branch estimator check, and a small real/complex layerwise gradient-access toy. These scripts make the construction inspectable without adding more detail to the manuscript. +- [`PAPER1_NUMERICS.md`](PAPER1_NUMERICS.md): guide to the first paper's numerical and reproducibility material. +- [`REPRODUCIBILITY.md`](REPRODUCIBILITY.md): commands for regenerating the first paper's datasets, diagnostics, safeguards, and figures. +- [`hopf_utils.py`](hopf_utils.py): shared real and complex Hopf maps, inverse maps, Jacobians, metrics, tangent-state assignments, and gate schedules. +- [`VQE_qibo.py`](VQE_qibo.py): the first paper's local real/complex VQE circuit-realizability demonstration. +- [`qbp_validation/`](qbp_validation/): exact-distribution Qibo validation for the second manuscript. -Hopf is used here as an optimization chart, not just as an amplitude loader. Existing state-preparation constructions address loading; this code focuses on the additional structures used in the paper's optimizer story: an explicit inverse chart, a diagonal metric, tangent-state assignments on the same circuit skeleton, and layerwise circuit checks. The finite-shot and Qibo/statevector demos are local estimator and circuit-realizability safeguards; they are not complete hardware-resource or total sampling-complexity analyses. - -## Synthetic tasks - -The main multi-size synthetic dataset uses six scrambled real-state tasks. Each task is scrambled by a fixed real orthogonal circuit, preventing the optimizers from exploiting an obvious computational-basis target. The focused complex experiment in `hopf_complex.py` repeats the same six objective families at `n=6` with complex unitary scrambling. - -The VQE test contains three scrambled Hamiltonian tests: - -**Parent Hamiltonian:** recover a scrambled target state by minimizing the parent-Hamiltonian gap. The task cost is the missing target overlap, - -```math -\mathcal{C}_{\mathrm{parent}}(\psi) -= -1 - |\langle \tau|\psi\rangle|^2 . -``` - -Here τ is the scrambled target state. The minimum is zero, reached when the variational state ψ equals τ up to the real global sign convention used in the code. - -**Hamming spectrum:** minimize a structured diagonal Hamming-distance spectrum hidden by a real orthogonal scrambling circuit. The task cost is the expected normalized Hamming distance in the unscrambled basis, - -```math -\mathcal{C}_{\mathrm{Ham}}(\psi) -= -\sum_x -\frac{d_{\mathrm{H}}(x,x_0)}{n} -\, -|(S^\top\psi)_x|^2 . -``` - -Here S is the fixed real scrambler and x₀ is the hidden target bit string. Physically, this tests whether the optimizer can find the scrambled computational-basis ground state of a simple but hidden diagonal spectrum. - -**Small-gap spectrum:** minimize a scrambled Hamiltonian with a deliberately small spectral gap near the ground state. The task cost is the expected value of a scrambled diagonal spectrum, - -```math -\mathcal{C}_{\mathrm{gap}}(\psi) -= -\sum_x -E_x -|(S^\top\psi)_x|^2 , -``` - -with - -```math -E_{x_0}=0,\qquad -E_{x_1}=10^{-2},\qquad -E_x=1\ \text{for all other }x . -``` - -The physical stress feature is the nearby distractor state x₁: the optimizer must distinguish the true ground state from a low-lying excited state. - -The metrology test contains three nonlinear sensing-inspired tests: - -**Single-target Fisher:** maximize a fixed-readout Fisher proxy for one scrambled target state. The task cost is - -```math -\mathcal{C}_{\mathrm{single}}(\psi) -= -1 - F, -\qquad -F = -|\langle \tau|\psi\rangle|^4 . -``` - -Equivalently, the code first computes the target probability A = |⟨τ|ψ⟩|² and then uses F = A². Physically, this rewards concentration of the probe state on one scrambled readout pattern. - -**QFI superposition:** maximize normalized pure-state QFI for a scrambled diagonal generator. In the unscrambled basis, the code computes - -```math -\mu = -\sum_x g_x |(S^\top\psi)_x|^2, -\qquad -\nu = -\sum_x g_x^2 |(S^\top\psi)_x|^2, -``` - -and then uses the normalized QFI objective - -```math -F_Q^{\mathrm{norm}}(\psi) -= -\frac{4(\nu-\mu^2)}{\mathrm{span}(G)^2}. -``` - -The task cost is - -```math -\mathcal{C}_{\mathrm{QFI}}(\psi) -= -1 - F_Q^{\mathrm{norm}}(\psi). -``` - -Physically, this rewards a probe state with large generator variance. The optimum is an equal superposition of the scrambled minimum- and maximum-generator eigenstates. - -**Balanced Fisher:** optimize a two-target Fisher objective where the optimum requires balanced overlap with two scrambled target states. The code computes - -```math -e_1 = |\langle \tau_1|\psi\rangle|^2, -\qquad -e_2 = |\langle \tau_2|\psi\rangle|^2, -``` - -then - -```math -F_1=e_1^2, -\qquad -F_2=e_2^2. -``` - -The balanced Fisher score is a soft minimum, - -```math -F_{\mathrm{bal}}(\psi) -= --\frac{1}{\beta} -\log\!\left( -\frac{ -e^{-\beta F_1} -+ -e^{-\beta F_2} -}{2} -\right), -\qquad -\beta=20 . -``` - -The task cost is - -```math -\mathcal{C}_{\mathrm{bal}}(\psi) -= -\frac{1}{4} -- -F_{\mathrm{bal}}(\psi). -``` - -The optimum has balanced probability on the two scrambled targets, giving approximately e₁ = e₂ = 1/2, hence F₁ = F₂ = 1/4 and zero cost. - -Lower cost or gap is better in all diagnostic summaries. - -To regenerate the real-Hopf plots locally, generate the datasets and run `plot_hopf.py`; see the sections below. The complex-Hopf script generates its CSV, diagnostics, and summary plot through one entry point. - -## Requirements +## Quick start Use Python 3.10 or newer. -Install the required packages for data generation, diagnostics, CNOT counting, and plotting: +### Paper I: numerical safeguards ```bash python -m pip install -r requirements.txt -``` - -Optional packages for explicit circuit-level checks: - -```bash python -m pip install -r requirements-optional.txt -``` - -`qibo` is used by optional circuit checks in `hopf_utils.py` and by `VQE_qibo.py` when running the explicit Qibo sampler. The synthetic data-generation, diagnostic, CNOT-count, finite-shot, and plotting scripts do not require Qibo. - -For a minimal local check without generating the full CSV archive, run the smoke-test commands in `REPRODUCIBILITY.md`. - -## Optimizer modes - -### Geometry-native Hopf optimizers - -Generated by `hopf_data.py`: - -```text -Hopf-EGT-CG -Hopf-Riemannian-LBFGS -Hopf-Riemannian-BB -``` - -These modes use the same cost-plus-Hopf-coordinate-gradient interface. The Hopf diagonal metric is used to lift coordinate gradients to state-sphere gradients, and accepted states are mapped back to Hopf coordinates. - -### Adam baselines - -Generated by `adam_data.py`: - -```text -Hopf-Adam -Mottonen-ideal-PS-Adam -``` - -`Hopf-Adam` applies Adam directly to Hopf coordinates. - -`Mottonen-ideal-PS-Adam` applies Adam to the physical post-multiplexing Möttönen rotation angles with an exact gradient equivalent to infinite-shot parameter shift. - -Both Adam baselines use adaptive cost-only backtracking. The diagnostics record the associated trial-evaluation counter, but these baselines are intended as implementation diagnostics rather than a standalone resource-accounting study. - -## Complex Hopf stress test - -`hopf_complex.py` provides a focused repository-level check of the complex Hopf extension. It repeats the same three VQE and three metrology-inspired objective families used in the real study, replacing the real orthogonal scrambling circuits with complex unitary scrambling circuits. It depends on `hopf_utils.py` in the same directory. - -The default experiment uses: - -```text -n = 6 -3 VQE tasks -3 metrology-inspired tasks -10 deterministic initial states per task -200 accepted optimizer updates -``` - -The four tracks are: - -```text -Hopf-Adam -Hopf-EGT-CG -Hopf-Riemannian-BB -Hopf-Riemannian-LBFGS -``` - -A complex Möttönen-Adam baseline is not included because a directly comparable extension of our chosen real post-multiplexing $R_y$ chart to a complex $R_y/R_z$ physical-angle chart was not implemented or validated here. - -Run the default experiment with: - -```bash -python hopf_complex.py -``` - -This writes: - -```text -complex_hopf_stress_data.csv -hopf_complex.png -``` - -The CSV is a generated artifact and is not intended to be tracked in source control. The PNG contains one two-panel final-gap figure with the same meaning, optimizer order, plotting conventions, and color code as the upper row of the main real-Hopf summary figure. - -For a minimal local check, run: - -```bash -python hopf_complex.py --quick -``` - -To regenerate the diagnostics and plot from an existing CSV without rerunning the optimizers: - -```bash -python hopf_complex.py --plot-only -``` - -The default run completed all `240` expected task-seed-mode traces and all `48,240` expected rows, with no missing steps, nonfinite gaps, materially negative gaps, or traces whose final gap was worse than their initial gap. The maximum state-norm error was `4.441e-16`. R-BB reached the `10^-8` threshold in all `30` VQE and all `30` metrology-inspired traces. EGT-CG and R-LBFGS had near-machine-precision medians with small slow-convergence tails, while coordinate Hopf-Adam retained substantially larger aggregate gaps. - -![View the complex-Hopf stress-test summary](hopf_complex.png) - -Roundoff-level negative final gaps are retained in the terminal diagnostics and clipped only when plotting on a logarithmic axis. - -## Seed convention - -By default, every fixed synthetic task is run from `10` deterministic initial-state seeds. - -The default seed rule is: - -```text -run_seed = problem_seed + seed_offset + seed_index -``` - -with: - -```text -seed_offset = 77777 -seed_index = 0, 1, ..., num_seeds - 1 -``` - -The default `num_seeds` is `10`. `hopf_complex.py` uses the same rule for its complex initial states. - -To change the number of initial states, pass the same value to both the data-generation and diagnostic scripts: - -```bash -python hopf_data.py --n 8 --num-seeds 5 --outdir data -python adam_data.py --n 8 --num-seeds 5 --outdir data - -python diagnose_hopf.py --indir data --ns 8 --num-seeds 5 -python diagnose_adam.py --indir data --ns 8 --num-seeds 5 -``` - -## Diagnostics - -Run diagnostics after generating the datasets. - -```bash -repo_root="$(pwd)" -mkdir -p diagnostics - -python diagnose_hopf.py \ - --indir data \ - --ns 6-10 \ - --steps 200 \ - --num-seeds 10 \ - --out "$repo_root/diagnostics/hopf_data_diagnostics.txt" - -python diagnose_adam.py \ - --indir data \ - --ns 6-10 \ - --steps 200 \ - --num-seeds 10 \ - --out "$repo_root/diagnostics/adam_data_diagnostics.txt" -``` - -The reports check: - -- expected files and row counts; -- expected task, mode, seed, and step completeness; -- CSV parse errors; -- sampled parameter and gradient vector lengths; -- nonfinite metrics or gradients; -- state-norm errors; -- aggregate final gaps; -- win counts and final rankings. - -The diagnostic scripts are streaming and avoid loading the full vector columns into pandas. -`diagnose_hopf.py` and `diagnose_adam.py` resolve relative `--out` paths inside `--indir`; the commands above use absolute output paths to keep the reports in the top-level `diagnostics/` directory. - -## Optional Hopf utility checks - -Run: - -```bash python hopf_utils.py +MPLBACKEND=Agg python VQE_qibo.py --steps 2 --shots 20 --sampler statevector --log-every 0 ``` -This executes consistency checks for inverse mapping, tangent-state synthesis, metric/Jacobian agreement, and optional Qibo circuit agreement when Qibo is installed. - -## Safeguard scripts - -The repository includes three standalone safeguard scripts. They are separate from the optimizer stress-test pipeline. `finite_shot_sanity.png` and `VQE_qibo.png` are included as small release panels for the estimator-level and circuit-level checks. These safeguards test specific parts of the construction; they are not full hardware-resource or total sampling-complexity benchmarks. - -### CNOT-count safeguard - -Run: - -```bash -python hopf_gate_count.py -``` - -This prints a table of real and complex Hopf CNOT counts for the default range `n = 4, ..., 20`, using the no-clean-ancilla CNOT model stated in the paper. It also saves a local plot named: - -```text -hopf_gate_count.pdf -``` - -To choose a smaller range and write the local output into a generated-figure folder: - -```bash -mkdir -p figures_safeguards - -python hopf_gate_count.py \ - --nmin 2 \ - --nmax 10 \ - --out figures_safeguards/hopf_gate_count.pdf -``` - -The script compares counts obtained directly from the generated Hopf gate schedules against the closed-form binomial count formulas. The saved plot is only a local check and is not shown in this README. +See [`PAPER1_NUMERICS.md`](PAPER1_NUMERICS.md) and +[`REPRODUCIBILITY.md`](REPRODUCIBILITY.md) for the full data-generation and plotting workflow. -### Finite-shot signed-branch gradient safeguard - -`finite_shot_sanity_check.py` isolates the statistical layer of the Hopf gradient protocol. It fixes one real Hopf state at `n=6`, uses a normalized diagonal Hamming-spectrum Hamiltonian, prepares the exact normalized tangent state for every one of the `63` Hopf coordinates, and estimates the symmetric signed-branch identity - -```math -\partial_{\theta_i} E -= -\sqrt{g_{i,i}} -\left(E_{\varphi_i^{(+)}}-E_{\varphi_i^{(-)}}\right). -``` - -The exact branch expectations first reproduce the independent tree-gradient reference to machine precision. Finite-shot branch energies are then sampled over `50` independent trials at each default shot count. - -Run: - -```bash -python finite_shot_sanity_check.py -``` - -This writes: - -```text -finite_shot_sanity_data.csv -finite_shot_sanity.png -``` - -To generate both PNG and PDF summaries, pass: +### Paper II: exact-logical Qibo validation ```bash -python finite_shot_sanity_check.py --formats png pdf -``` - -To regenerate the plot from the existing CSV without rerunning the sampling experiment: - -```bash -python finite_shot_sanity_check.py --plot-only -``` - -Default diagnostics: - -```text -exact branch-identity residual 2.513e-16 -Tree/Jacobian reference residual 2.355e-16 -maximum || 2.220e-16 - -shots/branch mean relative error SEM -100 0.186020 0.005300 -1,000 0.056339 0.001429 -10,000 0.018599 0.000477 - -log-log slope of mean error -0.500036 -ideal independent-shot slope -0.5 -``` - -The plotted error bars are standard errors over the `50` trials. The fitted slope is consistent with the expected `S^{-1/2}` statistical convergence. In the script's convention, `S` samples are allocated to each of the `+` and `-` branch-energy estimates for each component; this fixed-state diagonal-observable check does not model indexed-label allocation, Pauli grouping, device noise, or a global full-gradient shot budget. - -

- -

- - -### Qibo layerwise gradient-access safeguard - -Run: - -```bash -MPLBACKEND=Agg python VQE_qibo.py -``` - -This script runs two local `n=4` VQE toy models: - -- **Real Hopf:** a real nearest-neighbor chain with `X`, `Z`, `XX`, `YY`, and `ZZ` terms. -- **Complex Hopf:** a local chiral chain with `X`, `Y`, `Z`, `XX`, `YY`, `ZZ`, and `XY - YX` terms. - -For each toy, Adam uses a fixed learning rate and a fully random initialization. The plot compares the exact Hopf-gradient trajectory with the sampled layerwise-circuit trajectory. The layerwise estimator uses the baseline energy, indexed derivative-layer energies, and indexed signed-branch energies; the complex toy also uses the single indexed phase layer. The explicit-Qibo path builds label-controlled local circuits for clarity at `n=4`. Here `Ns` denotes the number of shots per indexed label/sign for each Hamiltonian Pauli term. - -The setting counts per Pauli readout are: - -```text -real Hopf: 1 + 2*4 = 9 -complex Hopf: 1 + 2*5 = 11 -``` - -The script writes: - -```text -VQE_qibo.png +python -m pip install -r requirements.txt +python -m pip install -r requirements-optional.txt +python -m unittest discover -s qbp_validation/tests -v ``` -![Qibo layerwise-circuit safeguard output](VQE_qibo.png) - -The included panel shows the sampled layerwise-circuit trajectories tracking the exact-gradient trajectories and approaching the exact ground-energy lines for both local toy Hamiltonians. This is a local circuit-realizability check for the Hopf gradient-access protocol, not a synthetic stress-test result or an asymptotic gate-count benchmark. - -## Reproducibility notes - -The scripts use deterministic random seeds by default. Re-running the same commands with the same Python, NumPy, and SciPy versions should reproduce the same synthetic task definitions and initial-state seeds. - -In Hopf CSVs, `last_line_evals` is a line-search work counter. It includes trial objective evaluations and, when strong-Wolfe curvature checks are reached, gradient-oracle calls used only by the line search. In Adam CSVs, `last_line_evals` counts adaptive cost-only trial evaluations. - -The multi-size datasets generated by `hopf_data.py` and `adam_data.py` use the real Hopf chart. `hopf_complex.py` provides a focused `n=6` complex-chart extension using the same six objective families and four Hopf optimizer tracks. `finite_shot_sanity_check.py` tests the finite-shot symmetric branch estimator at a fixed real state, while `VQE_qibo.py` exercises real and complex layerwise constructions in small local circuit-level examples. `hopf_utils.py` contains the shared real and complex Hopf maps, inverse maps, Jacobians, metrics, and tangent-state utilities. - -## Citation - -If you use this code, cite the accompanying Hopf ansatz paper. +The Paper-II suite evaluates exact Qibo statevectors and complete measurement +probability distributions for small systems. It creates no plots or benchmark datasets. +See [`qbp_validation/README.md`](qbp_validation/README.md) for the validation matrix and +scope of each test group. diff --git a/qbp_validation/README.md b/qbp_validation/README.md new file mode 100644 index 0000000..8544644 --- /dev/null +++ b/qbp_validation/README.md @@ -0,0 +1,104 @@ +# Exact-logical Qibo validation for Hopf quantum backpropagation + +This directory contains deterministic Qibo circuit-level validation for +*The Compass in Reverse: Quantum Backpropagation with the Hopf Ansatz*. +The manuscript is theoretical and makes no numerical-performance claim. The suite checks +exact-logical state, frame, measurement, and decoder identities using small statevector +circuits and independent analytic references. + +It covers the real global frame, complex magnitude and phase families, real and complex +magnitude checkpoints, singular chart cases, representative polyspherical trees, the +complete four-qubit construction, and the assigned compiler-ledger arithmetic. + +It does **not** benchmark optimization, finite-shot convergence, execution time, memory, +routing, approximate synthesis, device noise, or hardware performance. Compiler-ledger +arithmetic is checked separately and is not inferred from Qibo-transpiled gate counts. + +## What counts as circuit-level validation + +A circuit check in this directory satisfies both conditions: + +1. the tested state or measurement distribution is produced by a + `qibo.models.Circuit` executed with Qibo's NumPy statevector backend; +2. the expected result is computed independently by the recursive NumPy formulas in + `reference.py`, which contains no Qibo imports. + +The tests do not inject an expected statevector into Qibo. They build the Hopf frame and +breadth-first checkpoint circuits from controlled `RY(2 theta)` gates with explicit open +controls. The complex diagonal phase layer and phase-calibrated controlled observable are +represented as exact logical `Unitary` gates, matching the manuscript's access model. +Representative polyspherical pair rotations are also exact two-level `Unitary` gates; those +checks validate the finite-dimensional frame algebra, not an arbitrary-tree synthesis claim. + +All estimator checks use the full output statevector after the stated basis rotations and +sum the complete probability distribution exactly. No Monte Carlo shots are used. + +## Layout + +| File | Role | +|---|---| +| `conventions.py` | Big-endian basis labels, node/marker maps, gate specifications, polyspherical topology, and assigned ledger formulas. | +| `reference.py` | Independent recursive states, frame columns, metric factors, coordinate derivatives, and exact gradients. | +| `circuits.py` | Qibo builders for real frames, complex phase layers, global estimators, checkpoints, and polyspherical frames. | +| `decoders.py` | Walsh, signed-histogram, signed one-hot, and checkpoint decoders. | +| `cases.py` | Deterministic interior, singular, observable, and tree cases. | +| `tests/` | Circuit, decoder, singular-case, four-qubit, polyspherical, and ledger checks. | + +## Validation matrix + +| Test group | Manuscript statement checked | Qibo calculation | Independent reference | +|---|---|---|---| +| `test_conventions.py` | Basis ordering, marker labels, open controls, and the corrected `S^dagger`-then-`H` Y readout | Small deterministic circuits | Explicit bits and matrices | +| `test_real_frame.py` | `W_R |0> = |psi>` and `W_R |lambda(j)> = |e_j>` including singular completions | Qibo circuit unitary and statevector | Recursive state and local complements | +| `test_real_global_estimator.py` | One all-X distribution returns every real magnitude derivative | Exact Qibo probabilities | Analytic coordinate gradient | +| `test_complex_magnitude.py` | Phase-dressed frame columns and the asymmetric magnitude decoder | Qibo frame/phase circuits and probabilities | Complex magnitude derivatives | +| `test_complex_phase.py` | Ancilla-Y/system-Z signed one-hot phase estimator and common-phase cancellation | Exact Qibo probabilities | Analytic leaf-phase derivatives | +| `test_checkpoints.py` | Real and complex signed one-hot magnitude checkpoints at every depth | Qibo Y/Z/Y distributions | Analytic layer gradients | +| `test_singular_cases.py` | Zero metric factors give zero coordinate derivatives without breaking the frame | Boundary-angle circuits | Recursive zero derivatives | +| `test_polyspherical.py` | Local-complement frame, diagonal metric, translation, and parity estimator | Qibo two-level circuits | Independent ordered-tree recursion | +| `test_decoders.py` | FWHT, signed leaf histogram, suffix marginalization, and one-hot decoding | Deterministic probability arrays | Direct sums | +| `test_four_qubit_example.py` | All appendix markers, node-5 formulas, and all four checkpoint settings | Fixed `n=4` circuits | Displayed manuscript formulas | +| `test_resource_ledger.py` | Assigned controlled-rotation, frame, depth, and checkpoint charges | No transpiler count is used | Exact integer formulas | + +## Coverage + +The balanced Hopf checks use `n = 1, 2, 3, 4`. They include: + +- deterministic generic interior parameters; +- real final-layer angles beyond `pi/2`, exercising the sign convention; +- zero and `pi/2` upstream angles, producing zero metric factors; +- complex states with zero-amplitude leaves; +- Pauli, diagonal-reflection, and fixed-seed Householder observables satisfying + `O = O^dagger` and `O^2 = I`. + +The polyspherical checks use two ordered full binary trees on eight encoded leaves, including +an unbalanced tree with a nonzero root anchor so that the theorem's Pauli-X translation is +exercised explicitly. + +## Run + +From the repository root: + +```bash +python -m pip install -r requirements.txt +python -m pip install -r requirements-optional.txt +python -m unittest discover -s qbp_validation/tests -v +``` + +The suite was developed and checked with Qibo 0.3.4 and its NumPy backend. A separate pinned +environment, validation CLI, and GitHub Actions workflow are intentionally deferred. + +A successful run executes the complete deterministic test suite and creates no figures, +CSVs, or benchmark summaries. + +## Interpretation + +The tests validate finite-size instances of the exact algebraic premises used in the +manuscript. They do not numerically prove concentration inequalities or asymptotic resource +statements. In particular: + +- Hoeffding and fixed-norm concentration results remain mathematical theorems; +- decoder complexity and storage claims are not timed empirically; +- the polyspherical resource extension remains conditional on the stated efficient-frame + realization assumptions; +- Qibo gate counts are not compared with the manuscript's assigned compiler ledger. diff --git a/qbp_validation/__init__.py b/qbp_validation/__init__.py new file mode 100644 index 0000000..a4a2ba1 --- /dev/null +++ b/qbp_validation/__init__.py @@ -0,0 +1,3 @@ +"""Deterministic exact-logical Qibo validation for the Hopf-QBP manuscript.""" + +__all__ = ["cases", "circuits", "conventions", "decoders", "reference"] diff --git a/qbp_validation/cases.py b/qbp_validation/cases.py new file mode 100644 index 0000000..c416cad --- /dev/null +++ b/qbp_validation/cases.py @@ -0,0 +1,111 @@ +"""Deterministic validation cases and Hermitian-unitary observables.""" +from __future__ import annotations + +import numpy as np + +from .conventions import PolyBranch, PolyLeaf, PolyTree + + +def regular_magnitudes(n: int, seed: int = 1000) -> np.ndarray: + rng = np.random.default_rng(seed + n) + N = 1 << n + theta = rng.uniform(0.17, 1.31, size=N - 1) + # Exercise the real final-layer sign convention without approaching zeros. + final_start = (1 << (n - 1)) - 1 + theta[final_start:] = rng.uniform(0.2, 2.0 * np.pi - 0.2, size=N // 2) + return theta + + +def complex_magnitudes(n: int, seed: int = 2000) -> np.ndarray: + rng = np.random.default_rng(seed + n) + return rng.uniform(0.14, 1.40, size=(1 << n) - 1) + + +def phases(n: int, seed: int = 3000) -> np.ndarray: + rng = np.random.default_rng(seed + n) + return rng.uniform(0.0, 2.0 * np.pi, size=1 << n) + + +def singular_magnitudes(n: int) -> np.ndarray: + theta = complex_magnitudes(n, seed=4000) + if n == 1: + theta[0] = 0.0 + else: + theta[0] = 0.0 + theta[1] = np.pi / 2.0 + return theta + + +def pauli_matrix(letter: str) -> np.ndarray: + matrices = { + "I": np.eye(2, dtype=complex), + "X": np.asarray([[0, 1], [1, 0]], dtype=complex), + "Y": np.asarray([[0, -1j], [1j, 0]], dtype=complex), + "Z": np.asarray([[1, 0], [0, -1]], dtype=complex), + } + return matrices[letter] + + +def pauli_string(word: str) -> np.ndarray: + matrix = np.asarray([[1.0]], dtype=complex) + for letter in word: + matrix = np.kron(matrix, pauli_matrix(letter)) + return matrix + + +def diagonal_reflection(n: int) -> np.ndarray: + N = 1 << n + signs = np.asarray([1.0 if ((3 * index + index.bit_count()) % 5) < 3 else -1.0 for index in range(N)]) + return np.diag(signs.astype(complex)) + + +def householder_reflection(n: int, seed: int = 5000) -> np.ndarray: + rng = np.random.default_rng(seed + n) + N = 1 << n + vector = rng.normal(size=N) + 1j * rng.normal(size=N) + vector /= np.linalg.norm(vector) + return np.eye(N, dtype=complex) - 2.0 * np.outer(vector, np.conjugate(vector)) + + +def observables(n: int) -> tuple[np.ndarray, ...]: + if n == 1: + word = "Y" + else: + letters = "XYZI" + word = "".join(letters[index % len(letters)] for index in range(n)) + return (pauli_string(word), diagonal_reflection(n), householder_reflection(n)) + + +def balanced_poly_tree(labels: list[int], prefix: str = "v") -> PolyTree: + if len(labels) == 1: + return PolyLeaf(labels[0]) + midpoint = len(labels) // 2 + return PolyBranch( + f"{prefix}_{len(labels)}_{labels[0]}", + balanced_poly_tree(labels[:midpoint], prefix + "L"), + balanced_poly_tree(labels[midpoint:], prefix + "R"), + ) + + +def unbalanced_poly_tree() -> PolyTree: + # Eight leaves, deliberately permuted so the root anchor is nonzero and the + # Pauli-X translation in the theorem is exercised. + leaves = [PolyLeaf(label) for label in (3, 0, 6, 1, 7, 2, 5, 4)] + left = PolyBranch( + "u_left", + leaves[0], + PolyBranch("u_left_tail", leaves[1], PolyBranch("u_left_deep", leaves[2], leaves[3])), + ) + right = PolyBranch( + "u_right", + PolyBranch("u_right_a", leaves[4], leaves[5]), + PolyBranch("u_right_b", leaves[6], leaves[7]), + ) + return PolyBranch("u_root", left, right) + + +def poly_angles(tree: PolyTree, seed: int = 6000) -> dict[str, float]: + from .conventions import poly_preorder + + rng = np.random.default_rng(seed) + return {node.key: float(rng.uniform(-1.2, 1.2)) for node in poly_preorder(tree)} diff --git a/qbp_validation/circuits.py b/qbp_validation/circuits.py new file mode 100644 index 0000000..00d1c48 --- /dev/null +++ b/qbp_validation/circuits.py @@ -0,0 +1,351 @@ +"""Qibo circuit builders for deterministic exact-logical validation. + +The builders emit actual ``qibo.models.Circuit`` objects. Analytic state and +coordinate-derivative formulas live exclusively in :mod:`reference`. +""" +from __future__ import annotations + +from typing import Mapping, Sequence + +import numpy as np +from qibo import gates, models, set_backend + +from .conventions import ( + PolyBranch, + PolyTree, + bit_at, + depth_gate_specs, + frame_gate_specs, + infer_n_from_magnitudes, + poly_anchor, + poly_marker, + poly_preorder, + validate_poly_tree, +) + +set_backend("numpy") + +SDG = np.asarray([[1.0, 0.0], [0.0, -1.0j]], dtype=complex) + + +def _controlled_ry_with_pattern( + circuit: models.Circuit, + target: int, + controls: Sequence[int], + control_values: Sequence[int], + angle: float, +) -> None: + if len(controls) != len(control_values): + raise ValueError("Control qubits and values have different lengths.") + open_controls = [q for q, value in zip(controls, control_values) if value == 0] + for q in open_controls: + circuit.add(gates.X(q)) + gate = gates.RY(target, 2.0 * float(angle)) + if controls: + gate = gate.controlled_by(*controls) + circuit.add(gate) + for q in reversed(open_controls): + circuit.add(gates.X(q)) + + +def add_real_frame( + circuit: models.Circuit, + theta: np.ndarray, + system_qubits: Sequence[int], + *, + inverse: bool = False, +) -> None: + theta = np.asarray(theta, dtype=float).reshape(-1) + n = infer_n_from_magnitudes(theta) + if len(system_qubits) != n: + raise ValueError("System-qubit list length does not match theta.") + specs = list(frame_gate_specs(n)) + if inverse: + specs.reverse() + for spec in specs: + target = system_qubits[spec.target] + controls = [q for index, q in enumerate(system_qubits) if index != spec.target] + values = [bit_at(spec.anchor, index, n) for index in range(n) if index != spec.target] + angle = -theta[spec.node - 1] if inverse else theta[spec.node - 1] + _controlled_ry_with_pattern(circuit, target, controls, values, angle) + + +def add_depth_layer( + circuit: models.Circuit, + theta: np.ndarray, + system_qubits: Sequence[int], + depth: int, + *, + inverse: bool = False, +) -> None: + theta = np.asarray(theta, dtype=float).reshape(-1) + n = infer_n_from_magnitudes(theta) + if len(system_qubits) != n: + raise ValueError("System-qubit list length does not match theta.") + for spec in depth_gate_specs(n, depth): + controls = list(system_qubits[:depth]) + values = [((spec.position >> (depth - 1 - pos)) & 1) for pos in range(depth)] + angle = -theta[spec.node - 1] if inverse else theta[spec.node - 1] + _controlled_ry_with_pattern(circuit, system_qubits[depth], controls, values, angle) + + +def add_depth_preparation( + circuit: models.Circuit, + theta: np.ndarray, + system_qubits: Sequence[int], +) -> None: + n = infer_n_from_magnitudes(theta) + for depth in range(n): + add_depth_layer(circuit, theta, system_qubits, depth) + + +def add_inverse_depth_suffix( + circuit: models.Circuit, + theta: np.ndarray, + system_qubits: Sequence[int], + selected_depth: int, +) -> None: + n = infer_n_from_magnitudes(theta) + for depth in range(n - 1, selected_depth, -1): + add_depth_layer(circuit, theta, system_qubits, depth, inverse=True) + + +def add_phase_layer( + circuit: models.Circuit, + phases: np.ndarray, + system_qubits: Sequence[int], + *, + inverse: bool = False, +) -> None: + phases = np.asarray(phases, dtype=float).reshape(-1) + if phases.size != 1 << len(system_qubits): + raise ValueError("Phase vector length does not match the system register.") + sign = -1.0 if inverse else 1.0 + matrix = np.diag(np.exp(1j * sign * phases)) + circuit.add(gates.Unitary(matrix, *system_qubits)) + + +def add_controlled_observable( + circuit: models.Circuit, + observable: np.ndarray, + ancilla: int, + system_qubits: Sequence[int], +) -> None: + observable = np.asarray(observable, dtype=complex) + N = 1 << len(system_qubits) + if observable.shape != (N, N): + raise ValueError("Observable dimension does not match system register.") + circuit.add(gates.Unitary(observable, *system_qubits).controlled_by(ancilla)) + + +def add_x_basis_rotation(circuit: models.Circuit, qubits: Sequence[int]) -> None: + for q in qubits: + circuit.add(gates.H(q)) + + +def add_y_basis_rotation(circuit: models.Circuit, qubit: int) -> None: + # State-update order is S^dagger followed by H, i.e. net H S^dagger. + circuit.add(gates.Unitary(SDG, qubit)) + circuit.add(gates.H(qubit)) + + +def frame_circuit(theta: np.ndarray, *, inverse: bool = False) -> models.Circuit: + n = infer_n_from_magnitudes(theta) + circuit = models.Circuit(n) + add_real_frame(circuit, theta, tuple(range(n)), inverse=inverse) + return circuit + + +def depth_preparation_circuit(theta: np.ndarray) -> models.Circuit: + n = infer_n_from_magnitudes(theta) + circuit = models.Circuit(n) + add_depth_preparation(circuit, theta, tuple(range(n))) + return circuit + + +def real_global_measurement_circuit(theta: np.ndarray, observable: np.ndarray) -> models.Circuit: + n = infer_n_from_magnitudes(theta) + ancilla = 0 + system = tuple(range(1, n + 1)) + circuit = models.Circuit(n + 1) + circuit.add(gates.H(ancilla)) + add_real_frame(circuit, theta, system) + add_controlled_observable(circuit, observable, ancilla, system) + add_real_frame(circuit, theta, system, inverse=True) + add_x_basis_rotation(circuit, (ancilla, *system)) + return circuit + + +def complex_magnitude_measurement_circuit( + magnitudes: np.ndarray, + phases: np.ndarray, + observable: np.ndarray, +) -> models.Circuit: + n = infer_n_from_magnitudes(magnitudes) + ancilla = 0 + system = tuple(range(1, n + 1)) + circuit = models.Circuit(n + 1) + circuit.add(gates.H(ancilla)) + add_real_frame(circuit, magnitudes, system) + add_phase_layer(circuit, phases, system) + add_controlled_observable(circuit, observable, ancilla, system) + add_phase_layer(circuit, phases, system, inverse=True) + add_real_frame(circuit, magnitudes, system, inverse=True) + add_x_basis_rotation(circuit, (ancilla, *system)) + return circuit + + +def complex_phase_measurement_circuit( + magnitudes: np.ndarray, + phases: np.ndarray, + observable: np.ndarray, +) -> models.Circuit: + n = infer_n_from_magnitudes(magnitudes) + ancilla = 0 + system = tuple(range(1, n + 1)) + circuit = models.Circuit(n + 1) + circuit.add(gates.H(ancilla)) + add_real_frame(circuit, magnitudes, system) + add_phase_layer(circuit, phases, system) + add_controlled_observable(circuit, observable, ancilla, system) + add_y_basis_rotation(circuit, ancilla) + return circuit + + +def real_checkpoint_measurement_circuit( + theta: np.ndarray, + observable: np.ndarray, + depth: int, +) -> models.Circuit: + n = infer_n_from_magnitudes(theta) + if not 0 <= depth < n: + raise ValueError("Depth must lie in 0, ..., n-1.") + ancilla = 0 + system = tuple(range(1, n + 1)) + circuit = models.Circuit(n + 1) + circuit.add(gates.H(ancilla)) + add_depth_preparation(circuit, theta, system) + add_controlled_observable(circuit, observable, ancilla, system) + add_inverse_depth_suffix(circuit, theta, system, depth) + add_y_basis_rotation(circuit, ancilla) + add_y_basis_rotation(circuit, system[depth]) + return circuit + + +def complex_checkpoint_measurement_circuit( + magnitudes: np.ndarray, + phases: np.ndarray, + observable: np.ndarray, + depth: int, +) -> models.Circuit: + n = infer_n_from_magnitudes(magnitudes) + if not 0 <= depth < n: + raise ValueError("Depth must lie in 0, ..., n-1.") + ancilla = 0 + system = tuple(range(1, n + 1)) + circuit = models.Circuit(n + 1) + circuit.add(gates.H(ancilla)) + add_depth_preparation(circuit, magnitudes, system) + add_phase_layer(circuit, phases, system) + add_controlled_observable(circuit, observable, ancilla, system) + add_phase_layer(circuit, phases, system, inverse=True) + add_inverse_depth_suffix(circuit, magnitudes, system, depth) + add_y_basis_rotation(circuit, ancilla) + add_y_basis_rotation(circuit, system[depth]) + return circuit + + +def _two_level_rotation(dimension: int, first: int, second: int, angle: float) -> np.ndarray: + matrix = np.eye(dimension, dtype=complex) + c = float(np.cos(angle)) + s = float(np.sin(angle)) + matrix[first, first] = c + matrix[second, first] = s + matrix[first, second] = -s + matrix[second, second] = c + return matrix + + +def _add_poly_translation(circuit: models.Circuit, system_qubits: Sequence[int], label: int) -> None: + n = len(system_qubits) + for index, qubit in enumerate(system_qubits): + if bit_at(label, index, n): + circuit.add(gates.X(qubit)) + + +def add_polyspherical_frame( + circuit: models.Circuit, + tree: PolyTree, + angles: Mapping[str, float], + system_qubits: Sequence[int], + *, + shifted: bool = True, + inverse: bool = False, +) -> None: + n = len(system_qubits) + validate_poly_tree(tree, n, angles) + nodes = list(poly_preorder(tree)) + dimension = 1 << n + if not inverse: + if shifted: + _add_poly_translation(circuit, system_qubits, poly_anchor(tree)) + for node in nodes: + matrix = _two_level_rotation( + dimension, + poly_anchor(node.left), + poly_anchor(node.right), + float(angles[node.key]), + ) + circuit.add(gates.Unitary(matrix, *system_qubits)) + else: + for node in reversed(nodes): + matrix = _two_level_rotation( + dimension, + poly_anchor(node.left), + poly_anchor(node.right), + -float(angles[node.key]), + ) + circuit.add(gates.Unitary(matrix, *system_qubits)) + if shifted: + _add_poly_translation(circuit, system_qubits, poly_anchor(tree)) + + +def polyspherical_frame_circuit( + tree: PolyTree, + angles: Mapping[str, float], + n: int, + *, + shifted: bool = True, +) -> models.Circuit: + circuit = models.Circuit(n) + add_polyspherical_frame(circuit, tree, angles, tuple(range(n)), shifted=shifted) + return circuit + + +def polyspherical_global_measurement_circuit( + tree: PolyTree, + angles: Mapping[str, float], + n: int, + observable: np.ndarray, +) -> models.Circuit: + ancilla = 0 + system = tuple(range(1, n + 1)) + circuit = models.Circuit(n + 1) + circuit.add(gates.H(ancilla)) + add_polyspherical_frame(circuit, tree, angles, system, shifted=True) + add_controlled_observable(circuit, observable, ancilla, system) + add_polyspherical_frame(circuit, tree, angles, system, shifted=True, inverse=True) + add_x_basis_rotation(circuit, (ancilla, *system)) + return circuit + + +def statevector(circuit: models.Circuit) -> np.ndarray: + result = circuit() + return np.asarray(result.state(), dtype=complex).reshape(-1) + + +def probabilities(circuit: models.Circuit) -> np.ndarray: + state = statevector(circuit) + probs = np.abs(state) ** 2 + probs /= probs.sum() + return probs diff --git a/qbp_validation/conventions.py b/qbp_validation/conventions.py new file mode 100644 index 0000000..233764e --- /dev/null +++ b/qbp_validation/conventions.py @@ -0,0 +1,212 @@ +"""Shared indexing and tree conventions for the exact-logical validation suite. + +This module deliberately contains no Qibo imports. It fixes the manuscript's +big-endian basis labels and the combinatorial data used independently by the +reference and circuit implementations. +""" +from __future__ import annotations + +from dataclasses import dataclass +import math +from typing import Iterator, Mapping, TypeAlias + + +@dataclass(frozen=True) +class FrameGateSpec: + node: int + depth: int + position: int + anchor: int + marker: int + target: int # Qibo/system-list index, with 0 the most-significant qubit. + + +@dataclass(frozen=True) +class DepthGateSpec: + node: int + depth: int + position: int + target: int # Qibo/system-list index. + + +@dataclass(frozen=True) +class PolyLeaf: + """A leaf encoded by one computational-basis label.""" + + label: int + + +@dataclass(frozen=True) +class PolyBranch: + """An ordered full-binary-tree internal node. + + ``key`` identifies the angle in an external mapping. Keeping the topology + separate from the numerical angles makes it possible for the circuit and + reference paths to consume the same immutable test case without sharing + state-vector code. + """ + + key: str + left: "PolyTree" + right: "PolyTree" + + +PolyTree: TypeAlias = PolyLeaf | PolyBranch + + +def infer_n_from_magnitudes(theta: object) -> int: + length = len(theta) # type: ignore[arg-type] + n_float = math.log2(length + 1) + n = int(round(n_float)) + if n < 1 or (1 << n) - 1 != length: + raise ValueError("Magnitude vector must have length 2**n - 1.") + return n + + +def node_depth_position(node: int) -> tuple[int, int]: + if node < 1: + raise ValueError("Internal-node indices start at 1.") + depth = node.bit_length() - 1 + return depth, node - (1 << depth) + + +def anchor_label(node: int, n: int) -> int: + depth, position = node_depth_position(node) + if depth >= n: + raise ValueError("Node is not internal for the requested qubit count.") + return position << (n - depth) + + +def marker_label(node: int, n: int) -> int: + depth, position = node_depth_position(node) + if depth >= n: + raise ValueError("Node is not internal for the requested qubit count.") + return (2 * position + 1) << (n - depth - 1) + + +def frame_gate_specs(n: int) -> tuple[FrameGateSpec, ...]: + specs: list[FrameGateSpec] = [] + for depth in range(n): + for position in range(1 << depth): + node = (1 << depth) + position + specs.append( + FrameGateSpec( + node=node, + depth=depth, + position=position, + anchor=anchor_label(node, n), + marker=marker_label(node, n), + target=depth, + ) + ) + return tuple(specs) + + +def depth_gate_specs(n: int, depth: int | None = None) -> tuple[DepthGateSpec, ...]: + depths = range(n) if depth is None else (depth,) + specs: list[DepthGateSpec] = [] + for d in depths: + if not 0 <= d < n: + raise ValueError("Depth must lie in 0, ..., n-1.") + for position in range(1 << d): + specs.append( + DepthGateSpec( + node=(1 << d) + position, + depth=d, + position=position, + target=d, + ) + ) + return tuple(specs) + + +def bit_at(label: int, qubit: int, n: int) -> int: + """Return the bit on big-endian qubit ``qubit`` (0 is most significant).""" + + if not 0 <= qubit < n: + raise ValueError("Qubit index out of range.") + return (label >> (n - 1 - qubit)) & 1 + + +def parity(label: int, outcome: int) -> int: + """Return the mod-two inner product of two basis labels.""" + + return (label & outcome).bit_count() & 1 + + +def marker_map(n: int) -> dict[int, int]: + return {spec.node: spec.marker for spec in frame_gate_specs(n)} + + +def controlled_ry_cnot_charge(num_controls: int) -> int: + if num_controls < 0: + raise ValueError("Number of controls must be nonnegative.") + if num_controls == 0: + return 0 + if num_controls <= 4: + return (1 << (num_controls + 1)) - 2 + return 16 * (num_controls + 1) - 40 + + +def frame_cnot_charge(n: int) -> int: + return ((1 << n) - 1) * controlled_ry_cnot_charge(n - 1) + + +def depth_layer_cnot_charge(depth: int) -> int: + if depth < 0: + raise ValueError("Depth must be nonnegative.") + return (1 << depth) * controlled_ry_cnot_charge(depth) + + +def depth_preparation_cnot_charge(n: int) -> int: + return sum(depth_layer_cnot_charge(d) for d in range(n)) + + +def checkpoint_cnot_charge_without_observable(n: int, depth: int) -> int: + if not 0 <= depth < n: + raise ValueError("Depth must lie in 0, ..., n-1.") + suffix = sum(depth_layer_cnot_charge(d) for d in range(depth + 1, n)) + return depth_preparation_cnot_charge(n) + suffix + + +def poly_preorder(tree: PolyTree) -> Iterator[PolyBranch]: + if isinstance(tree, PolyLeaf): + return + yield tree + yield from poly_preorder(tree.left) + yield from poly_preorder(tree.right) + + +def poly_leaves(tree: PolyTree) -> tuple[int, ...]: + if isinstance(tree, PolyLeaf): + return (tree.label,) + return poly_leaves(tree.left) + poly_leaves(tree.right) + + +def poly_anchor(tree: PolyTree) -> int: + if isinstance(tree, PolyLeaf): + return tree.label + return poly_anchor(tree.left) + + +def poly_marker(node: PolyBranch) -> int: + return poly_anchor(node.right) + + +def poly_relative_markers(tree: PolyTree) -> dict[str, int]: + root_anchor = poly_anchor(tree) + return {node.key: root_anchor ^ poly_marker(node) for node in poly_preorder(tree)} + + +def validate_poly_tree(tree: PolyTree, n: int, angles: Mapping[str, float]) -> None: + labels = poly_leaves(tree) + if len(labels) != len(set(labels)): + raise ValueError("Polyspherical leaf labels must be distinct.") + if any(label < 0 or label >= (1 << n) for label in labels): + raise ValueError("Polyspherical leaf label outside the encoded space.") + keys = [node.key for node in poly_preorder(tree)] + if len(keys) != len(set(keys)): + raise ValueError("Polyspherical internal-node keys must be unique.") + missing = set(keys).difference(angles) + if missing: + raise ValueError(f"Missing angles for internal nodes: {sorted(missing)}") diff --git a/qbp_validation/decoders.py b/qbp_validation/decoders.py new file mode 100644 index 0000000..825f86c --- /dev/null +++ b/qbp_validation/decoders.py @@ -0,0 +1,115 @@ +"""Exact-distribution and record-level decoders for Hopf-QBP validation.""" +from __future__ import annotations + +import numpy as np + +from .conventions import marker_label, parity + + +def fwht(values: np.ndarray) -> np.ndarray: + """Unnormalized in-place-style Walsh-Hadamard transform on a copy.""" + + out = np.asarray(values, dtype=float).reshape(-1).copy() + n = out.size + if n == 0 or n & (n - 1): + raise ValueError("FWHT input length must be a positive power of two.") + width = 1 + while width < n: + for start in range(0, n, 2 * width): + a = out[start : start + width].copy() + b = out[start + width : start + 2 * width].copy() + out[start : start + width] = a + b + out[start + width : start + 2 * width] = a - b + width *= 2 + return out + + +def reshape_ancilla_system(probabilities: np.ndarray) -> np.ndarray: + probs = np.asarray(probabilities, dtype=float).reshape(-1) + if probs.size < 4 or probs.size & (probs.size - 1): + raise ValueError("Probability vector must describe one ancilla and at least one system qubit.") + N = probs.size // 2 + return probs.reshape(2, N) + + +def signed_system_histogram(probabilities: np.ndarray) -> np.ndarray: + probs = reshape_ancilla_system(probabilities) + return probs[0] - probs[1] + + +def global_moments_direct(probabilities: np.ndarray) -> np.ndarray: + signed = signed_system_histogram(probabilities) + N = signed.size + moments = np.zeros(N, dtype=float) + for label in range(N): + moments[label] = sum( + value * (-1.0 if parity(label, outcome) else 1.0) + for outcome, value in enumerate(signed) + ) + return moments + + +def global_moments_fwht(probabilities: np.ndarray) -> np.ndarray: + return fwht(signed_system_histogram(probabilities)) + + +def decode_balanced_magnitude_gradient( + probabilities: np.ndarray, + metric: np.ndarray, + n: int, + *, + use_fwht: bool = True, +) -> np.ndarray: + metric = np.asarray(metric, dtype=float).reshape(-1) + if metric.size != (1 << n) - 1: + raise ValueError("Metric length does not match n.") + moments = global_moments_fwht(probabilities) if use_fwht else global_moments_direct(probabilities) + return np.asarray( + [2.0 * np.sqrt(max(metric[node - 1], 0.0)) * moments[marker_label(node, n)] for node in range(1, 1 << n)], + dtype=float, + ) + + +def decode_phase_gradient(probabilities: np.ndarray) -> np.ndarray: + probs = reshape_ancilla_system(probabilities) + return 2.0 * (probs[0] - probs[1]) + + +def decode_checkpoint_gradient(probabilities: np.ndarray, n: int, depth: int) -> np.ndarray: + if not 0 <= depth < n: + raise ValueError("Depth must lie in 0, ..., n-1.") + probs = reshape_ancilla_system(probabilities) + width = 1 << depth + result = np.zeros(width, dtype=float) + target_shift = n - depth - 1 + for ancilla in (0, 1): + for system_label in range(1 << n): + prefix = system_label >> (n - depth) if depth else 0 + target = (system_label >> target_shift) & 1 + sign = -2.0 * (-1.0 if (ancilla + target) & 1 else 1.0) + result[prefix] += sign * probs[ancilla, system_label] + return result + + +def decode_polyspherical_gradient( + probabilities: np.ndarray, + incoming: dict[str, float], + relative_markers: dict[str, int], +) -> dict[str, float]: + moments = global_moments_fwht(probabilities) + return { + key: 2.0 * incoming[key] * moments[relative_markers[key]] + for key in relative_markers + } + + +def phase_record(ancilla_bit: int, leaf: int, N: int) -> np.ndarray: + record = np.zeros(N, dtype=float) + record[leaf] = 2.0 * (-1.0 if ancilla_bit else 1.0) + return record + + +def checkpoint_record(ancilla_bit: int, target_bit: int, prefix: int, width: int) -> np.ndarray: + record = np.zeros(width, dtype=float) + record[prefix] = -2.0 * (-1.0 if (ancilla_bit + target_bit) & 1 else 1.0) + return record diff --git a/qbp_validation/reference.py b/qbp_validation/reference.py new file mode 100644 index 0000000..b418733 --- /dev/null +++ b/qbp_validation/reference.py @@ -0,0 +1,280 @@ +"""Independent NumPy reference formulas for the Hopf-QBP manuscript. + +No function in this module imports Qibo or calls a circuit builder. The test +suite uses these recursive formulas as the analytic side of every circuit-level +comparison. +""" +from __future__ import annotations + +from dataclasses import dataclass +from typing import Mapping + +import numpy as np + +from .conventions import ( + PolyBranch, + PolyLeaf, + PolyTree, + infer_n_from_magnitudes, + marker_label, + poly_anchor, + poly_leaves, + poly_marker, + poly_preorder, + validate_poly_tree, +) + + +@dataclass(frozen=True) +class RealTreeData: + n: int + state: np.ndarray + subtree: tuple[np.ndarray | None, ...] + complements: tuple[np.ndarray | None, ...] + incoming: np.ndarray + metric: np.ndarray + derivatives: tuple[np.ndarray, ...] + + +@dataclass(frozen=True) +class PolyTreeData: + state: np.ndarray + complements: dict[str, np.ndarray] + incoming: dict[str, float] + derivatives: dict[str, np.ndarray] + anchors: dict[str, int] + markers: dict[str, int] + + +def basis_vector(dimension: int, label: int, *, dtype: type = complex) -> np.ndarray: + vec = np.zeros(dimension, dtype=dtype) + vec[label] = 1 + return vec + + +def real_tree_data(theta: np.ndarray | list[float] | tuple[float, ...]) -> RealTreeData: + theta = np.asarray(theta, dtype=float).reshape(-1) + n = infer_n_from_magnitudes(theta) + N = 1 << n + + subtree: list[np.ndarray | None] = [None] * (2 * N) + complements: list[np.ndarray | None] = [None] * N + for ell in range(N): + subtree[N + ell] = basis_vector(N, ell, dtype=float) + + for node in range(N - 1, 0, -1): + left = subtree[2 * node] + right = subtree[2 * node + 1] + assert left is not None and right is not None + c = float(np.cos(theta[node - 1])) + s = float(np.sin(theta[node - 1])) + subtree[node] = c * left + s * right + complements[node] = -s * left + c * right + + incoming = np.zeros(N, dtype=float) + incoming[1] = 1.0 + for node in range(1, N): + left = 2 * node + right = left + 1 + if left < N: + incoming[left] = incoming[node] * np.cos(theta[node - 1]) + if right < N: + incoming[right] = incoming[node] * np.sin(theta[node - 1]) + + derivatives: list[np.ndarray] = [] + for node in range(1, N): + comp = complements[node] + assert comp is not None + derivatives.append(incoming[node] * comp) + + root = subtree[1] + assert root is not None + return RealTreeData( + n=n, + state=np.asarray(root, dtype=float), + subtree=tuple(subtree), + complements=tuple(complements), + incoming=incoming, + metric=incoming[1:] ** 2, + derivatives=tuple(derivatives), + ) + + +def real_state(theta: np.ndarray | list[float] | tuple[float, ...]) -> np.ndarray: + return real_tree_data(theta).state + + +def real_frame_matrix(theta: np.ndarray | list[float] | tuple[float, ...]) -> np.ndarray: + data = real_tree_data(theta) + N = 1 << data.n + frame = np.zeros((N, N), dtype=float) + frame[:, 0] = data.state + for node in range(1, N): + comp = data.complements[node] + assert comp is not None + frame[:, marker_label(node, data.n)] = comp + return frame + + +def complex_state(magnitudes: np.ndarray, phases: np.ndarray) -> np.ndarray: + magnitudes = np.asarray(magnitudes, dtype=float).reshape(-1) + phases = np.asarray(phases, dtype=float).reshape(-1) + data = real_tree_data(magnitudes) + N = 1 << data.n + if phases.size != N: + raise ValueError("Complex Hopf chart requires one phase per leaf.") + return np.exp(1j * phases) * data.state + + +def complex_magnitude_derivatives(magnitudes: np.ndarray, phases: np.ndarray) -> tuple[np.ndarray, ...]: + data = real_tree_data(magnitudes) + phase = np.exp(1j * np.asarray(phases, dtype=float).reshape(-1)) + if phase.size != (1 << data.n): + raise ValueError("Complex Hopf chart requires one phase per leaf.") + return tuple(phase * deriv for deriv in data.derivatives) + + +def complex_phase_derivatives(magnitudes: np.ndarray, phases: np.ndarray) -> tuple[np.ndarray, ...]: + psi = complex_state(magnitudes, phases) + N = psi.size + out: list[np.ndarray] = [] + for ell in range(N): + deriv = np.zeros(N, dtype=complex) + deriv[ell] = 1j * psi[ell] + out.append(deriv) + return tuple(out) + + +def expectation(state: np.ndarray, observable: np.ndarray) -> float: + return float(np.real(np.vdot(state, observable @ state))) + + +def coordinate_gradient(state: np.ndarray, derivatives: tuple[np.ndarray, ...], observable: np.ndarray) -> np.ndarray: + Opsi = observable @ state + return np.asarray([2.0 * np.real(np.vdot(deriv, Opsi)) for deriv in derivatives], dtype=float) + + +def real_gradient(theta: np.ndarray, observable: np.ndarray) -> np.ndarray: + data = real_tree_data(theta) + return coordinate_gradient(data.state.astype(complex), data.derivatives, observable) + + +def complex_magnitude_gradient(magnitudes: np.ndarray, phases: np.ndarray, observable: np.ndarray) -> np.ndarray: + psi = complex_state(magnitudes, phases) + return coordinate_gradient(psi, complex_magnitude_derivatives(magnitudes, phases), observable) + + +def complex_phase_gradient(magnitudes: np.ndarray, phases: np.ndarray, observable: np.ndarray) -> np.ndarray: + psi = complex_state(magnitudes, phases) + return coordinate_gradient(psi, complex_phase_derivatives(magnitudes, phases), observable) + + +def complex_full_gradient(magnitudes: np.ndarray, phases: np.ndarray, observable: np.ndarray) -> np.ndarray: + return np.concatenate( + [ + complex_magnitude_gradient(magnitudes, phases, observable), + complex_phase_gradient(magnitudes, phases, observable), + ] + ) + + +def is_hermitian_unitary(observable: np.ndarray, atol: float = 1e-12) -> bool: + observable = np.asarray(observable, dtype=complex) + eye = np.eye(observable.shape[0], dtype=complex) + return bool( + observable.ndim == 2 + and observable.shape[0] == observable.shape[1] + and np.allclose(observable, observable.conj().T, atol=atol, rtol=0.0) + and np.allclose(observable @ observable, eye, atol=atol, rtol=0.0) + ) + + +def _poly_reference_node( + tree: PolyTree, + angles: Mapping[str, float], + dimension: int, + incoming: float, + data: dict[str, object], +) -> np.ndarray: + if isinstance(tree, PolyLeaf): + return basis_vector(dimension, tree.label, dtype=float) + + angle = float(angles[tree.key]) + left = _poly_reference_node(tree.left, angles, dimension, incoming * np.cos(angle), data) + right = _poly_reference_node(tree.right, angles, dimension, incoming * np.sin(angle), data) + c = float(np.cos(angle)) + s = float(np.sin(angle)) + state = c * left + s * right + comp = -s * left + c * right + + incoming_map = data["incoming"] + complements = data["complements"] + derivatives = data["derivatives"] + anchors = data["anchors"] + markers = data["markers"] + assert isinstance(incoming_map, dict) + assert isinstance(complements, dict) + assert isinstance(derivatives, dict) + assert isinstance(anchors, dict) + assert isinstance(markers, dict) + incoming_map[tree.key] = float(incoming) + complements[tree.key] = comp + derivatives[tree.key] = float(incoming) * comp + anchors[tree.key] = poly_anchor(tree) + markers[tree.key] = poly_marker(tree) + return state + + +def polyspherical_tree_data(tree: PolyTree, angles: Mapping[str, float], n: int) -> PolyTreeData: + validate_poly_tree(tree, n, angles) + dimension = 1 << n + data: dict[str, object] = { + "incoming": {}, + "complements": {}, + "derivatives": {}, + "anchors": {}, + "markers": {}, + } + state = _poly_reference_node(tree, angles, dimension, 1.0, data) + return PolyTreeData( + state=state, + complements=data["complements"], # type: ignore[arg-type] + incoming=data["incoming"], # type: ignore[arg-type] + derivatives=data["derivatives"], # type: ignore[arg-type] + anchors=data["anchors"], # type: ignore[arg-type] + markers=data["markers"], # type: ignore[arg-type] + ) + + +def polyspherical_frame_matrix(tree: PolyTree, angles: Mapping[str, float], n: int) -> np.ndarray: + data = polyspherical_tree_data(tree, angles, n) + dimension = 1 << n + root_anchor = poly_anchor(tree) + frame = np.eye(dimension, dtype=float) + encoded_labels = set(poly_leaves(tree)) + # Columns outside the encoded leaf set remain the computational basis. + for label in encoded_labels: + frame[:, label] = 0.0 + frame[:, root_anchor] = data.state + for node in poly_preorder(tree): + frame[:, poly_marker(node)] = data.complements[node.key] + return frame + + +def polyspherical_shifted_frame_matrix(tree: PolyTree, angles: Mapping[str, float], n: int) -> np.ndarray: + frame = polyspherical_frame_matrix(tree, angles, n) + dimension = 1 << n + anchor = poly_anchor(tree) + translation = np.zeros((dimension, dimension), dtype=float) + for label in range(dimension): + translation[label ^ anchor, label] = 1.0 + return frame @ translation + + +def polyspherical_gradient(tree: PolyTree, angles: Mapping[str, float], n: int, observable: np.ndarray) -> dict[str, float]: + data = polyspherical_tree_data(tree, angles, n) + Opsi = observable @ data.state + return { + node.key: float(2.0 * np.real(np.vdot(data.derivatives[node.key], Opsi))) + for node in poly_preorder(tree) + } diff --git a/qbp_validation/tests/__init__.py b/qbp_validation/tests/__init__.py new file mode 100644 index 0000000..493d36a --- /dev/null +++ b/qbp_validation/tests/__init__.py @@ -0,0 +1 @@ +"""Unit tests for the exact-logical Hopf-QBP validation suite.""" diff --git a/qbp_validation/tests/test_checkpoints.py b/qbp_validation/tests/test_checkpoints.py new file mode 100644 index 0000000..7a2a896 --- /dev/null +++ b/qbp_validation/tests/test_checkpoints.py @@ -0,0 +1,59 @@ +from __future__ import annotations + +import unittest + +import numpy as np + +from qbp_validation.cases import complex_magnitudes, observables, phases, regular_magnitudes +from qbp_validation.circuits import ( + complex_checkpoint_measurement_circuit, + probabilities, + real_checkpoint_measurement_circuit, +) +from qbp_validation.decoders import checkpoint_record, decode_checkpoint_gradient +from qbp_validation.reference import complex_magnitude_gradient, real_gradient + + +class CheckpointTests(unittest.TestCase): + def test_real_checkpoint_all_depths(self) -> None: + for n in range(1, 5): + theta = regular_magnitudes(n) + for obs_index, observable in enumerate(observables(n)): + exact = real_gradient(theta, observable) + for depth in range(n): + with self.subTest(n=n, observable=obs_index, depth=depth): + probs = probabilities(real_checkpoint_measurement_circuit(theta, observable, depth)) + decoded = decode_checkpoint_gradient(probs, n, depth) + start = (1 << depth) - 1 + stop = (1 << (depth + 1)) - 1 + np.testing.assert_allclose(decoded, exact[start:stop], atol=3e-12, rtol=0.0) + + def test_complex_magnitude_checkpoint_all_depths(self) -> None: + for n in range(1, 5): + magnitudes = complex_magnitudes(n) + leaf_phases = phases(n) + for obs_index, observable in enumerate(observables(n)): + exact = complex_magnitude_gradient(magnitudes, leaf_phases, observable) + for depth in range(n): + with self.subTest(n=n, observable=obs_index, depth=depth): + probs = probabilities( + complex_checkpoint_measurement_circuit( + magnitudes, leaf_phases, observable, depth + ) + ) + decoded = decode_checkpoint_gradient(probs, n, depth) + start = (1 << depth) - 1 + stop = (1 << (depth + 1)) - 1 + np.testing.assert_allclose(decoded, exact[start:stop], atol=3e-12, rtol=0.0) + + def test_checkpoint_records_have_fixed_norm_two(self) -> None: + for width in (1, 2, 4, 8): + for ancilla in (0, 1): + for target in (0, 1): + for prefix in range(width): + record = checkpoint_record(ancilla, target, prefix, width) + self.assertAlmostEqual(float(np.linalg.norm(record)), 2.0) + + +if __name__ == "__main__": + unittest.main() diff --git a/qbp_validation/tests/test_complex_magnitude.py b/qbp_validation/tests/test_complex_magnitude.py new file mode 100644 index 0000000..7eca550 --- /dev/null +++ b/qbp_validation/tests/test_complex_magnitude.py @@ -0,0 +1,60 @@ +from __future__ import annotations + +import unittest + +import numpy as np +from qibo import models + +from qbp_validation.cases import complex_magnitudes, observables, phases +from qbp_validation.circuits import ( + add_phase_layer, + add_real_frame, + complex_magnitude_measurement_circuit, + probabilities, +) +from qbp_validation.conventions import marker_label +from qbp_validation.decoders import decode_balanced_magnitude_gradient +from qbp_validation.reference import ( + complex_magnitude_derivatives, + complex_magnitude_gradient, + complex_state, + real_tree_data, +) + + +class ComplexMagnitudeTests(unittest.TestCase): + def test_phase_dressed_frame_columns(self) -> None: + for n in range(1, 5): + magnitudes = complex_magnitudes(n) + leaf_phases = phases(n) + circuit = models.Circuit(n) + system = tuple(range(n)) + add_real_frame(circuit, magnitudes, system) + add_phase_layer(circuit, leaf_phases, system) + unitary = np.asarray(circuit.unitary()) + np.testing.assert_allclose(unitary[:, 0], complex_state(magnitudes, leaf_phases), atol=2e-12, rtol=0.0) + derivatives = complex_magnitude_derivatives(magnitudes, leaf_phases) + incoming = real_tree_data(magnitudes).incoming + for node, derivative in enumerate(derivatives, start=1): + expected_column = derivative / incoming[node] + np.testing.assert_allclose( + unitary[:, marker_label(node, n)], expected_column, atol=2e-12, rtol=0.0 + ) + + def test_exact_magnitude_estimator(self) -> None: + for n in range(1, 5): + magnitudes = complex_magnitudes(n) + leaf_phases = phases(n) + metric = real_tree_data(magnitudes).metric + for obs_index, observable in enumerate(observables(n)): + with self.subTest(n=n, observable=obs_index): + probs = probabilities( + complex_magnitude_measurement_circuit(magnitudes, leaf_phases, observable) + ) + decoded = decode_balanced_magnitude_gradient(probs, metric, n) + exact = complex_magnitude_gradient(magnitudes, leaf_phases, observable) + np.testing.assert_allclose(decoded, exact, atol=3e-12, rtol=0.0) + + +if __name__ == "__main__": + unittest.main() diff --git a/qbp_validation/tests/test_complex_phase.py b/qbp_validation/tests/test_complex_phase.py new file mode 100644 index 0000000..0bc9c59 --- /dev/null +++ b/qbp_validation/tests/test_complex_phase.py @@ -0,0 +1,47 @@ +from __future__ import annotations + +import unittest + +import numpy as np + +from qbp_validation.cases import complex_magnitudes, observables, phases +from qbp_validation.circuits import complex_phase_measurement_circuit, probabilities +from qbp_validation.decoders import decode_phase_gradient, phase_record +from qbp_validation.reference import complex_phase_gradient + + +class ComplexPhaseTests(unittest.TestCase): + def test_signed_one_hot_phase_estimator(self) -> None: + for n in range(1, 5): + magnitudes = complex_magnitudes(n) + leaf_phases = phases(n) + for obs_index, observable in enumerate(observables(n)): + with self.subTest(n=n, observable=obs_index): + probs = probabilities(complex_phase_measurement_circuit(magnitudes, leaf_phases, observable)) + decoded = decode_phase_gradient(probs) + exact = complex_phase_gradient(magnitudes, leaf_phases, observable) + np.testing.assert_allclose(decoded, exact, atol=3e-12, rtol=0.0) + self.assertAlmostEqual(float(decoded.sum()), 0.0, places=11) + + def test_every_phase_record_has_fixed_norm_two(self) -> None: + for N in (2, 4, 8, 16): + for ancilla in (0, 1): + for leaf in range(N): + self.assertAlmostEqual(float(np.linalg.norm(phase_record(ancilla, leaf, N))), 2.0) + + def test_zero_amplitude_leaf_phase_derivatives_vanish(self) -> None: + n = 3 + magnitudes = complex_magnitudes(n) + magnitudes[0] = 0.0 # The entire left half has nonzero mass; right half vanishes. + leaf_phases = phases(n) + observable = observables(n)[-1] + decoded = decode_phase_gradient( + probabilities(complex_phase_measurement_circuit(magnitudes, leaf_phases, observable)) + ) + exact = complex_phase_gradient(magnitudes, leaf_phases, observable) + np.testing.assert_allclose(decoded, exact, atol=3e-12, rtol=0.0) + np.testing.assert_allclose(decoded[4:], 0.0, atol=3e-12, rtol=0.0) + + +if __name__ == "__main__": + unittest.main() diff --git a/qbp_validation/tests/test_conventions.py b/qbp_validation/tests/test_conventions.py new file mode 100644 index 0000000..083b148 --- /dev/null +++ b/qbp_validation/tests/test_conventions.py @@ -0,0 +1,54 @@ +from __future__ import annotations + +import unittest + +import numpy as np +from qibo import gates, models + +from qbp_validation.cases import observables +from qbp_validation.circuits import add_y_basis_rotation, probabilities +from qbp_validation.conventions import bit_at, frame_gate_specs, marker_label, parity +from qbp_validation.reference import is_hermitian_unitary + + +class ConventionTests(unittest.TestCase): + def test_marker_map_is_bijection_onto_nonzero_labels(self) -> None: + for n in range(1, 6): + labels = [marker_label(node, n) for node in range(1, 1 << n)] + self.assertEqual(sorted(labels), list(range(1, 1 << n))) + + def test_big_endian_bit_convention(self) -> None: + label = 0b1011 + self.assertEqual([bit_at(label, q, 4) for q in range(4)], [1, 0, 1, 1]) + self.assertEqual(parity(0b0110, 0b1011), 1) + + def test_qibo_qubit_zero_is_most_significant(self) -> None: + circuit = models.Circuit(3) + circuit.add(gates.X(0)) + state = np.asarray(circuit().state()) + self.assertEqual(int(np.argmax(np.abs(state))), 0b100) + + def test_y_basis_rotation_has_manuscript_sign(self) -> None: + # Prepare |+i> = S H |0>. S^dagger followed by H must map it to |0>. + circuit = models.Circuit(1) + circuit.add(gates.H(0)) + circuit.add(gates.Unitary(np.diag([1.0, 1.0j]), 0)) + add_y_basis_rotation(circuit, 0) + probs = probabilities(circuit) + np.testing.assert_allclose(probs, [1.0, 0.0], atol=1e-13, rtol=0.0) + + def test_frame_specs_match_active_pairs(self) -> None: + for n in range(1, 5): + for spec in frame_gate_specs(n): + self.assertEqual(spec.anchor ^ spec.marker, 1 << (n - 1 - spec.target)) + self.assertEqual(bit_at(spec.anchor, spec.target, n), 0) + self.assertEqual(bit_at(spec.marker, spec.target, n), 1) + + def test_validation_observables_satisfy_access_model(self) -> None: + for n in range(1, 5): + for observable in observables(n): + self.assertTrue(is_hermitian_unitary(observable)) + + +if __name__ == "__main__": + unittest.main() diff --git a/qbp_validation/tests/test_decoders.py b/qbp_validation/tests/test_decoders.py new file mode 100644 index 0000000..1e75955 --- /dev/null +++ b/qbp_validation/tests/test_decoders.py @@ -0,0 +1,57 @@ +from __future__ import annotations + +import unittest + +import numpy as np + +from qbp_validation.decoders import ( + decode_checkpoint_gradient, + decode_phase_gradient, + fwht, + global_moments_direct, + global_moments_fwht, +) + + +class DecoderTests(unittest.TestCase): + def test_fwht_matches_dense_walsh_matrix(self) -> None: + for n in range(1, 7): + N = 1 << n + rng = np.random.default_rng(7000 + n) + values = rng.normal(size=N) + dense = np.asarray( + [ + [(-1.0) ** ((row & column).bit_count() & 1) for column in range(N)] + for row in range(N) + ] + ) + np.testing.assert_allclose(fwht(values), dense @ values, atol=2e-12, rtol=0.0) + + def test_global_direct_and_fwht_decoders_agree(self) -> None: + for n in range(1, 6): + rng = np.random.default_rng(7100 + n) + probs = rng.uniform(size=2 * (1 << n)) + probs /= probs.sum() + np.testing.assert_allclose( + global_moments_direct(probs), global_moments_fwht(probs), atol=2e-12, rtol=0.0 + ) + + def test_phase_decoder_is_signed_leaf_histogram(self) -> None: + probs = np.asarray([[0.10, 0.20, 0.05, 0.15], [0.05, 0.10, 0.25, 0.10]]) + np.testing.assert_allclose(decode_phase_gradient(probs.reshape(-1)), 2.0 * (probs[0] - probs[1])) + + def test_checkpoint_decoder_discards_suffix_bits(self) -> None: + n = 3 + depth = 1 + probs = np.zeros((2, 1 << n), dtype=float) + # Same ancilla/prefix/target data but two different suffix values. + probs[0, 0b010] = 0.2 + probs[0, 0b011] = 0.3 + probs[1, 0b110] = 0.1 + probs[1, 0b111] = 0.4 + decoded = decode_checkpoint_gradient(probs.reshape(-1), n, depth) + np.testing.assert_allclose(decoded, [1.0, -1.0], atol=1e-13, rtol=0.0) + + +if __name__ == "__main__": + unittest.main() diff --git a/qbp_validation/tests/test_four_qubit_example.py b/qbp_validation/tests/test_four_qubit_example.py new file mode 100644 index 0000000..9770261 --- /dev/null +++ b/qbp_validation/tests/test_four_qubit_example.py @@ -0,0 +1,60 @@ +from __future__ import annotations + +import unittest + +import numpy as np + +from qbp_validation.cases import observables, regular_magnitudes +from qbp_validation.circuits import probabilities, real_checkpoint_measurement_circuit +from qbp_validation.conventions import marker_map +from qbp_validation.decoders import decode_checkpoint_gradient +from qbp_validation.reference import real_gradient, real_tree_data + + +class FourQubitExampleTests(unittest.TestCase): + def test_all_appendix_marker_labels(self) -> None: + expected = { + 1: 0b1000, + 2: 0b0100, + 3: 0b1100, + 4: 0b0010, + 5: 0b0110, + 6: 0b1010, + 7: 0b1110, + 8: 0b0001, + 9: 0b0011, + 10: 0b0101, + 11: 0b0111, + 12: 0b1001, + 13: 0b1011, + 14: 0b1101, + 15: 0b1111, + } + self.assertEqual(marker_map(4), expected) + + def test_explicit_node_five_column_and_metric(self) -> None: + theta = regular_magnitudes(4) + data = real_tree_data(theta) + expected = np.zeros(16, dtype=float) + expected[0b0100] = -np.sin(theta[4]) * np.cos(theta[9]) + expected[0b0101] = -np.sin(theta[4]) * np.sin(theta[9]) + expected[0b0110] = np.cos(theta[4]) * np.cos(theta[10]) + expected[0b0111] = np.cos(theta[4]) * np.sin(theta[10]) + np.testing.assert_allclose(data.complements[5], expected, atol=2e-12, rtol=0.0) + self.assertAlmostEqual(data.incoming[5], np.cos(theta[0]) * np.sin(theta[1]), places=12) + + def test_all_four_checkpoint_settings(self) -> None: + theta = regular_magnitudes(4) + observable = observables(4)[-1] + exact = real_gradient(theta, observable) + for depth in range(4): + decoded = decode_checkpoint_gradient( + probabilities(real_checkpoint_measurement_circuit(theta, observable, depth)), 4, depth + ) + start = (1 << depth) - 1 + stop = (1 << (depth + 1)) - 1 + np.testing.assert_allclose(decoded, exact[start:stop], atol=3e-12, rtol=0.0) + + +if __name__ == "__main__": + unittest.main() diff --git a/qbp_validation/tests/test_polyspherical.py b/qbp_validation/tests/test_polyspherical.py new file mode 100644 index 0000000..b9267b4 --- /dev/null +++ b/qbp_validation/tests/test_polyspherical.py @@ -0,0 +1,77 @@ +from __future__ import annotations + +import unittest + +import numpy as np + +from qbp_validation.cases import ( + balanced_poly_tree, + observables, + poly_angles, + unbalanced_poly_tree, +) +from qbp_validation.circuits import ( + polyspherical_frame_circuit, + polyspherical_global_measurement_circuit, + probabilities, +) +from qbp_validation.conventions import poly_relative_markers +from qbp_validation.decoders import decode_polyspherical_gradient +from qbp_validation.reference import ( + polyspherical_gradient, + polyspherical_shifted_frame_matrix, + polyspherical_tree_data, +) + + +class PolysphericalTests(unittest.TestCase): + def _cases(self): + return ( + (balanced_poly_tree(list(range(8))), poly_angles(balanced_poly_tree(list(range(8))), 6100)), + (unbalanced_poly_tree(), poly_angles(unbalanced_poly_tree(), 6200)), + ) + + def test_frame_and_local_complements(self) -> None: + n = 3 + for case_index, (tree, angles) in enumerate(self._cases()): + with self.subTest(case=case_index): + qibo_frame = np.asarray(polyspherical_frame_circuit(tree, angles, n).unitary()) + reference = polyspherical_shifted_frame_matrix(tree, angles, n) + np.testing.assert_allclose(qibo_frame, reference, atol=3e-12, rtol=0.0) + np.testing.assert_allclose( + qibo_frame.conj().T @ qibo_frame, np.eye(1 << n), atol=3e-12, rtol=0.0 + ) + + def test_parity_estimator_for_representative_trees(self) -> None: + n = 3 + for case_index, (tree, angles) in enumerate(self._cases()): + data = polyspherical_tree_data(tree, angles, n) + for obs_index, observable in enumerate(observables(n)): + with self.subTest(case=case_index, observable=obs_index): + probs = probabilities( + polyspherical_global_measurement_circuit(tree, angles, n, observable) + ) + decoded = decode_polyspherical_gradient( + probs, data.incoming, poly_relative_markers(tree) + ) + exact = polyspherical_gradient(tree, angles, n, observable) + for key in exact: + self.assertAlmostEqual(decoded[key], exact[key], places=11) + + def test_derivative_factorization_and_diagonal_metric(self) -> None: + tree = unbalanced_poly_tree() + angles = poly_angles(tree, 6300) + data = polyspherical_tree_data(tree, angles, 3) + keys = list(data.derivatives) + gram = np.asarray( + [ + [np.real(np.vdot(data.derivatives[a], data.derivatives[b])) for b in keys] + for a in keys + ] + ) + expected = np.diag([data.incoming[key] ** 2 for key in keys]) + np.testing.assert_allclose(gram, expected, atol=3e-12, rtol=0.0) + + +if __name__ == "__main__": + unittest.main() diff --git a/qbp_validation/tests/test_real_frame.py b/qbp_validation/tests/test_real_frame.py new file mode 100644 index 0000000..71cc4e4 --- /dev/null +++ b/qbp_validation/tests/test_real_frame.py @@ -0,0 +1,50 @@ +from __future__ import annotations + +import unittest + +import numpy as np + +from qbp_validation.cases import regular_magnitudes, singular_magnitudes +from qbp_validation.circuits import depth_preparation_circuit, frame_circuit, statevector +from qbp_validation.conventions import marker_label +from qbp_validation.reference import real_frame_matrix, real_tree_data + + +class RealFrameTests(unittest.TestCase): + def _check_case(self, theta: np.ndarray) -> None: + data = real_tree_data(theta) + N = 1 << data.n + qibo_unitary = np.asarray(frame_circuit(theta).unitary()) + reference = real_frame_matrix(theta) + np.testing.assert_allclose(qibo_unitary, reference, atol=2e-12, rtol=0.0) + np.testing.assert_allclose(qibo_unitary.conj().T @ qibo_unitary, np.eye(N), atol=2e-12, rtol=0.0) + np.testing.assert_allclose(qibo_unitary[:, 0], data.state, atol=2e-12, rtol=0.0) + + for node in range(1, N): + marker = marker_label(node, data.n) + complement = data.complements[node] + self.assertIsNotNone(complement) + np.testing.assert_allclose(qibo_unitary[:, marker], complement, atol=2e-12, rtol=0.0) + np.testing.assert_allclose( + data.derivatives[node - 1], + data.incoming[node] * complement, + atol=2e-12, + rtol=0.0, + ) + + depth_state = statevector(depth_preparation_circuit(theta)) + np.testing.assert_allclose(depth_state, data.state, atol=2e-12, rtol=0.0) + + def test_regular_frames_n1_to_n4(self) -> None: + for n in range(1, 5): + with self.subTest(n=n): + self._check_case(regular_magnitudes(n)) + + def test_singular_frame_completion_remains_unitary(self) -> None: + for n in range(1, 5): + with self.subTest(n=n): + self._check_case(singular_magnitudes(n)) + + +if __name__ == "__main__": + unittest.main() diff --git a/qbp_validation/tests/test_real_global_estimator.py b/qbp_validation/tests/test_real_global_estimator.py new file mode 100644 index 0000000..313205f --- /dev/null +++ b/qbp_validation/tests/test_real_global_estimator.py @@ -0,0 +1,43 @@ +from __future__ import annotations + +import unittest + +import numpy as np + +from qbp_validation.cases import observables, regular_magnitudes +from qbp_validation.circuits import probabilities, real_global_measurement_circuit +from qbp_validation.decoders import ( + decode_balanced_magnitude_gradient, + global_moments_direct, + global_moments_fwht, +) +from qbp_validation.reference import real_gradient, real_tree_data + + +class RealGlobalEstimatorTests(unittest.TestCase): + def test_exact_distribution_returns_complete_gradient(self) -> None: + for n in range(1, 5): + theta = regular_magnitudes(n) + data = real_tree_data(theta) + for obs_index, observable in enumerate(observables(n)): + with self.subTest(n=n, observable=obs_index): + probs = probabilities(real_global_measurement_circuit(theta, observable)) + self.assertAlmostEqual(float(probs.sum()), 1.0, places=13) + np.testing.assert_allclose( + global_moments_fwht(probs), + global_moments_direct(probs), + atol=2e-12, + rtol=0.0, + ) + decoded = decode_balanced_magnitude_gradient(probs, data.metric, n) + exact = real_gradient(theta, observable) + np.testing.assert_allclose(decoded, exact, atol=3e-12, rtol=0.0) + + def test_per_record_coordinate_range_is_at_most_two(self) -> None: + for n in range(1, 6): + metric = real_tree_data(regular_magnitudes(n)).metric + self.assertLessEqual(float(np.max(2.0 * np.sqrt(metric))), 2.0 + 1e-13) + + +if __name__ == "__main__": + unittest.main() diff --git a/qbp_validation/tests/test_resource_ledger.py b/qbp_validation/tests/test_resource_ledger.py new file mode 100644 index 0000000..f99cd1d --- /dev/null +++ b/qbp_validation/tests/test_resource_ledger.py @@ -0,0 +1,32 @@ +from __future__ import annotations + +import unittest + +from qbp_validation.conventions import ( + checkpoint_cnot_charge_without_observable, + controlled_ry_cnot_charge, + depth_layer_cnot_charge, + depth_preparation_cnot_charge, + frame_cnot_charge, +) + + +class ResourceLedgerTests(unittest.TestCase): + def test_assigned_controlled_ry_charges(self) -> None: + expected = {0: 0, 1: 2, 2: 6, 3: 14, 4: 30, 5: 56, 6: 72, 7: 88} + self.assertEqual({m: controlled_ry_cnot_charge(m) for m in expected}, expected) + + def test_four_qubit_frame_and_depth_totals(self) -> None: + self.assertEqual(frame_cnot_charge(4), 210) + self.assertEqual([depth_layer_cnot_charge(d) for d in range(4)], [0, 4, 24, 112]) + self.assertEqual(depth_preparation_cnot_charge(4), 140) + + def test_four_qubit_checkpoint_totals(self) -> None: + self.assertEqual( + [checkpoint_cnot_charge_without_observable(4, d) for d in range(4)], + [280, 276, 252, 140], + ) + + +if __name__ == "__main__": + unittest.main() diff --git a/qbp_validation/tests/test_singular_cases.py b/qbp_validation/tests/test_singular_cases.py new file mode 100644 index 0000000..ac4ca79 --- /dev/null +++ b/qbp_validation/tests/test_singular_cases.py @@ -0,0 +1,60 @@ +from __future__ import annotations + +import unittest + +import numpy as np + +from qbp_validation.cases import observables, phases, singular_magnitudes +from qbp_validation.circuits import ( + complex_checkpoint_measurement_circuit, + complex_magnitude_measurement_circuit, + probabilities, + real_checkpoint_measurement_circuit, + real_global_measurement_circuit, +) +from qbp_validation.decoders import decode_balanced_magnitude_gradient, decode_checkpoint_gradient +from qbp_validation.reference import complex_magnitude_gradient, real_gradient, real_tree_data + + +class SingularChartTests(unittest.TestCase): + def test_zero_metric_coordinates_decode_to_zero(self) -> None: + for n in range(2, 5): + theta = singular_magnitudes(n) + data = real_tree_data(theta) + zero = np.flatnonzero(data.metric < 1e-28) + self.assertGreater(zero.size, 0) + observable = observables(n)[-1] + probs = probabilities(real_global_measurement_circuit(theta, observable)) + decoded = decode_balanced_magnitude_gradient(probs, data.metric, n) + exact = real_gradient(theta, observable) + np.testing.assert_allclose(decoded, exact, atol=3e-12, rtol=0.0) + np.testing.assert_allclose(decoded[zero], 0.0, atol=3e-12, rtol=0.0) + + def test_real_and_complex_checkpoints_at_singular_points(self) -> None: + for n in range(2, 5): + theta = singular_magnitudes(n) + leaf_phases = phases(n) + observable = observables(n)[1] + exact_real = real_gradient(theta, observable) + exact_complex = complex_magnitude_gradient(theta, leaf_phases, observable) + for depth in range(n): + start = (1 << depth) - 1 + stop = (1 << (depth + 1)) - 1 + real_decoded = decode_checkpoint_gradient( + probabilities(real_checkpoint_measurement_circuit(theta, observable, depth)), n, depth + ) + complex_decoded = decode_checkpoint_gradient( + probabilities( + complex_checkpoint_measurement_circuit(theta, leaf_phases, observable, depth) + ), + n, + depth, + ) + np.testing.assert_allclose(real_decoded, exact_real[start:stop], atol=3e-12, rtol=0.0) + np.testing.assert_allclose( + complex_decoded, exact_complex[start:stop], atol=3e-12, rtol=0.0 + ) + + +if __name__ == "__main__": + unittest.main()