diff --git a/benchmarks/README.md b/benchmarks/README.md index 39dc22a30..65b95b269 100644 --- a/benchmarks/README.md +++ b/benchmarks/README.md @@ -87,3 +87,19 @@ uv run --with matplotlib python benchmarks/plot_branch_coalesce.py \ - `plot_branch_coalesce.py` — left panel: time vs `m` (log-log); right panel: sort-merge speedup `t_hash / t_sortmerge` vs `m`, with the crossover line and the "hash wins" band. + +# Cross-library Pauli propagation + +`ppvm` against [PauliPropagation.jl][xpp], [PauliStrings.jl][xps], Qiskit's +[pauli-prop][xqk] and Algorithmiq's [monoprop][xmp], on TFIM Trotter and +Heisenberg autocorrelator workloads, all single-threaded and all validated to +propagate the identical operator term-for-term before anything is timed. + +Lives in [`cross-library/`](cross-library/README.md) — see that README for the +circuit spec, the shared environment contract, the per-engine caveats (notably +that monoprop is parallel unless capped), and a measured run. + +[xpp]: https://github.com/MSRudolph/PauliPropagation.jl +[xps]: https://github.com/nicolasloizeau/PauliStrings.jl +[xqk]: https://github.com/Qiskit/pauli-prop +[xmp]: https://github.com/Algorithmiq/monoprop diff --git a/benchmarks/cross-library/README.md b/benchmarks/cross-library/README.md new file mode 100644 index 000000000..4c02f6232 --- /dev/null +++ b/benchmarks/cross-library/README.md @@ -0,0 +1,382 @@ +# Cross-library Pauli-propagation benchmark + +`ppvm` against the four other single-threaded Pauli-propagation engines we know +of, on two workloads: + +| library | language | entry point used | +|---|---|---| +| **`ppvm`** (this repo) | Rust | `PauliSum>` | +| [PauliPropagation.jl][pp] | Julia | `propagate(PauliRotation(...), psum; min_abs_coeff)` | +| [PauliStrings.jl][ps] | Julia | `trotter_step!(O, gates; truncation, truncate_every)` | +| [pauli-prop][qk] (Qiskit) | Rust-accelerated Python | `propagate_through_circuit(op, qc, max_terms, atol, frame="h")` | +| [monoprop][mp] (Algorithmiq) | C++ with Python bindings | `PauliPropagator.from_circuit(circuit, op, cutoff, lower_atol)` | + +```bash +# Everything, with the term-for-term agreement check first. +uv run --no-project python3 benchmarks/cross-library/run_xbench.py \ + --qubits-tfim 8,16,24,32,40,48,56,64 \ + --qubits-heisenberg 6,8,10,12,14 \ + --steps 10 --atol 1e-6 --iters 2 --out target/xbench + +uv run --no-project --with matplotlib python3 benchmarks/cross-library/plot_xbench.py \ + --csv target/xbench/results.csv --out target/xbench/xbench.png \ + --title "Pauli propagation: ppvm vs PauliPropagation.jl, PauliStrings.jl, pauli-prop, monoprop" + +# How wrong each engine is at a given atol, rather than how fast — see +# "Known differences" below. Regenerates accuracy.csv. +uv run --no-project python3 benchmarks/cross-library/xbench_accuracy.py \ + --out benchmarks/cross-library/accuracy.csv + +# The regime where the two truncation rules actually diverge. Three sweeps make +# up accuracy_divergence.csv — depth, angle, and a scrambling circuit with small +# angles and great depth (plus a θ=π/4 control on the same family). +A=benchmarks/cross-library/xbench_accuracy.py +uv run --no-project python3 $A --models heisenberg --libs ppvm,monoprop \ + --qubits 6 --dt 0.05 --steps 10,20,40,80,160,320 --atols 1e-3,1e-5 --out d.csv +uv run --no-project python3 $A --models heisenberg --libs ppvm,monoprop \ + --qubits 8 --steps 10 --dt 0.02,0.05,0.1,0.2,0.4,0.6,0.785 --atols 1e-3 --out a.csv +uv run --no-project python3 $A --models scramble --libs ppvm,monoprop --qubits 8 \ + --steps 400 --dt 0.05 --seeds 1,2,3,4,5 --atols 1e-3,1e-4,1e-5 --out s.csv + +# Four panels: atol scaling, the depth trend, the angle reversal, the scrambler. +uv run --no-project --with matplotlib python3 benchmarks/cross-library/plot_accuracy.py \ + --csv benchmarks/cross-library/accuracy.csv \ + --divergence benchmarks/cross-library/accuracy_divergence.csv \ + --out target/xbench/accuracy.png +``` + +A subset, if you prefer: `--libs ppvm,monoprop`. + +## What is measured + +Both workloads are **Heisenberg-picture propagation of an observable through an +explicit first-order Trotter product of Pauli rotations**, with a +coefficient-magnitude truncation after every gate. Everything below is fixed +across all five engines — the gate list, its order, the angles, the truncation +rule, and the readout. + +### `tfim` — transverse-field Ising, magnetization + +`H = J Σᵢ ZᵢZᵢ₊₁ + h Σᵢ Xᵢ` on an open chain. One step, in this order: + +1. `RX(2h·dt)` on site `0, 1, …, n−1` +2. `RZZ(2J·dt)` on bond `(0,1), (1,2), …, (n−2,n−1)` + +Observable `O = Σᵢ Zᵢ`; readout `⟨0…0|O(t)|0…0⟩`, i.e. the sum of the +coefficients of the X-free terms. + +### `heisenberg` — isotropic Heisenberg + field, autocorrelator + +`H = J Σᵢ (XᵢXᵢ₊₁ + YᵢYᵢ₊₁ + ZᵢZᵢ₊₁) + h Σᵢ Zᵢ`. One step, in this order: + +1. `RXX`, `RYY`, `RZZ` at `2J·dt` on bond `(0,1)`, then `(1,2)`, … +2. `RZ(2h·dt)` on site `0, 1, …, n−1` + +Observable `O = Z₀`; readout the autocorrelator `S(t) = tr[Z₀·O(t)]/2ⁿ`, which +is just the coefficient of `Z₀` (the Paulis are orthonormal under that pairing). + +### `scramble` — random all-to-all rotations, autocorrelator + +`steps·n` two-qubit Pauli rotations `exp(−iθ/2·Pₐ⊗P_b)`, each drawing a uniformly +random **all-to-all** pair `a≠b`, random axes `Pₐ, P_b ∈ {X,Y,Z}` and a random +angle `θ ∈ (0, 2J·dt]`. Observable and readout as for `heisenberg`. + +Only `ppvm` and `monoprop` implement it (`--libs ppvm,monoprop`); it exists +because the two Trotter models are a weak test bed for *truncation* questions. +They are nearest-neighbour, uniform-angle and conserve enough structure that +their support saturates a symmetry sector — at `n=8` Heisenberg reaches exactly +16 384 of 65 536 words, and the other 49 152 are zero by symmetry rather than by +dynamics. `scramble` has no lattice, no conserved quantity and no repeated angle: +at `n=8` it fills all 65 535 non-identity words, and since unitary conjugation +preserves the Frobenius norm the whole coefficient vector has `‖c‖₂ = 1` spread +across them. That makes it the right instance for asking whether a truncation +rule can resolve one coefficient against a scrambled background. + +Both runners generate the circuit from a **splitmix64 stream reimplemented +identically in Rust and Python**, seeded by `SEED`. Neither language's stdlib RNG +is specified tightly enough to rely on, and the draw order (pair, offset, axis, +axis, angle) has to match exactly — including the `b = (a+1+r mod n−1) mod n` +trick, which avoids a rejection loop that would consume a variable number of +draws and desynchronise the two streams. The term-for-term dump diff is what +enforces this; it catches any divergence immediately. + +`θ = 2·c·dt` for a Hamiltonian term `c·G` is the convention four of the five +engines use directly for `exp(iθ/2·G)·P·exp(−iθ/2·G)`, so the propagated +operators are identical, not merely similar. monoprop is the exception and needs +a conversion; see below. + +## The parameter contract + +Every runner reads the same environment variables and writes the same CSV, so +they can be run directly as well as through the driver: + +| variable | meaning | +|---|---| +| `MODEL` | `tfim`, `heisenberg`, or `scramble` (`ppvm` and `monoprop` only) | +| `QUBITS` | comma-separated widths | +| `STEPS` | Trotter steps | +| `DT`, `JCOUP`, `HFIELD` | `dt`, `J`, `h` | +| `ATOL` | truncation threshold on `|c|` | +| `ITERS` | timed repeats; the **minimum** is reported | +| `SEED` | `scramble` only — the circuit seed | +| `DUMP` | print the propagated support instead of timing it | +| `MAX_TERMS` | `pauli-prop` only — its mandatory cap (see below) | +| `monoprop_NUM_THREADS`, `monoprop_PARTITIONS` | `monoprop` only — the thread cap (see below) | + +CSV columns: `model,library,qubits,steps,dt,atol,time_s,terms,observable`. + +```bash +MODEL=heisenberg QUBITS=8,10,12 STEPS=10 ATOL=1e-6 ITERS=2 \ + cargo run --release -p ppvm-pauli-sum --example xbench +``` + +## Validation — why the driver refuses to time first + +`run_xbench.py` runs every engine with `DUMP=1` at `n=4, steps=3, atol=1e-14` +and diffs the whole propagated support against `ppvm`'s, term for term. It +aborts on any missing term, extra term, or coefficient difference above +`--validate-tol` (1e-10 by default). + +This is not ceremony. Three real bugs in this harness produced numbers that +looked completely reasonable and were wrong: + +* **The `pauli-prop` circuit was reversed.** In the Heisenberg frame it + conjugates from the *end* of the instruction list backwards, so the spec's + gate order needs `reversed()` on append. Appending forward propagates a + different operator — 108 terms against TFIM's 124, 61 against Heisenberg's 64, + coefficients off by up to 0.1 — while the readout still matched to 9 digits. +* **Duplicate Paulis were silently collapsing.** `propagate_through_circuit` + returns a `SparsePauliOp` that may list the same Pauli more than once, so + `len(op)` is a row count, not a support size, and a `{label: coeff}` dict + comprehension keeps the last duplicate instead of summing them. +* **monoprop's angle convention is neither the spec's nor Qiskit's.** Its + `ExpGate` applies `exp(+iθH)`, so the spec's `exp(−iθ_spec/2·G)` needs + `θ = −θ_spec/2`, *and* the gate list needs reversing like `pauli-prop`'s. All + eight sign/order combinations were tried against the reference: forward order + loses terms outright (109 against TFIM's 124), and every wrong angle keeps the + right support while moving coefficients by up to 1.3. Only `θ = −θ_spec/2` on + a reversed list lands within 5e-13. + +An observable is one scalar and can agree by luck or by cancellation. A +term-for-term diff cannot. + +## monoprop is parallel unless you stop it + +This is the one caveat that changes a headline number rather than a decimal, so +it gets its own section. + +monoprop takes **one serial partition per physical core** when left alone: with a +single MPI rank, `resolve_partition_count_` reads the core count and fans out. +The PyPI wheels are built without MPI (`monoprop.has_mpi == False`), which is +easy to misread as "therefore serial" — it is not. On the 14-core machine below, +Heisenberg at `n=12` measures: + +| | wall | CPU (user+sys) | CPU/wall | +|---|---:|---:|---:| +| default | 0.66 s | 6.28 s | 9.48× | +| `monoprop_NUM_THREADS=1 monoprop_PARTITIONS=off` | 1.97 s | 1.97 s | 1.00× | + +So an uncapped monoprop reports a **3× faster** wall time than the serial one, +against engines that never had the option. Both variables are read once into a +cached C++ static, so they must be in the environment before the first +propagation — the runner sets them itself as well as receiving them from the +driver. + +Because a stale environment variable fails silently, `xbench_monoprop.py` +**measures its own CPU/wall ratio around the timed region and exits non-zero** if +it exceeds `CPU_WALL_MAX` (1.5). Every monoprop row in the run below reported +`cpu/wall=1.00`. + +## Known differences between the engines + +These are real and are the reason the plot shows the **workload size** next to +the runtime. Read them before quoting a ratio. + +* **`pauli-prop` and `monoprop` truncate by a different rule.** Both prune at + branch-creation time, where the other three accumulate first and drop the + merged coefficient. `ppvm`, PauliPropagation.jl and PauliStrings.jl then agree + on the support **exactly**, at every width, on both models; the other two do + not. On the sweep below `pauli-prop` runs 11–36 % *under* the reference support + on TFIM (the gap widening with `n`) and 4–15 % *over* it on Heisenberg; + `monoprop` is within 1 % on TFIM and 4–12 % over on Heisenberg. The driver + prints a `note:` line for any engine more than 2 % off, and the plot's third + panel is there so the ratio is never read without it. +* **Prune-at-creation costs a little accuracy, not an order of it.** Measured as + the `L2` distance of the whole coefficient vector from a converged reference + (`ppvm` at `atol=1e-16`; the Heisenberg `n=8` sector saturates at 16 384 terms, + so that reference is exact), across `atol` from `1e-3` to `1e-7` at `n=8`, + `steps=10`: `monoprop` lands within 1.00–1.12× of `ppvm`'s error and + `pauli-prop` at 1.23–2.60×. Error falls by a decade for every decade of `atol` + in all five engines, so the two policies share a convergence order and differ + only in the constant. Don't read this off the `observable` column instead — it + is one scalar and its truncation error changes sign, so `ppvm` on TFIM goes + `7.7e-4, 1.8e-3, 1.0e-4, 3.0e-7, 1.2e-6` over those five thresholds. That dip + is a cancellation, and comparing against it manufactures a 27× gap where the + vector norm shows 6 %. `xbench_accuracy.py` writes both, for that reason; + [`accuracy.csv`](accuracy.csv) is the run these figures come from. +* **`monoprop` tracks more rows than it reports.** It retains monomials whose + coefficient has cancelled to exactly zero. This is the visible end of + prune-at-creation: a branch that cleared the threshold when it was emitted is + never re-tested after later contributions cancel it. At Heisenberg `n=8`, + `atol=1e-4` it holds 2 074 terms whose converged value is below `atol`, against + `ppvm`'s 1 107 — so roughly 7 % of its support is dead weight. That is a cost + in work rather than in accuracy; the coefficients it keeps are still right. Its `size()` at Heisenberg `n=14` is + 7 355 928 rows against the 3 204 697 terms above threshold — 2.3× — while on + TFIM the two are within 1 %. The `terms` column is the above-threshold support, + which is what the other four engines mean by it; `size()` goes to stderr next + to each row. +* **`pauli-prop` also has a mandatory `max_terms` cap** with pre-allocation. The + runner defaults it to `2²²` and **fails the run** if the support ever reaches + it, rather than quietly reporting a differently-truncated number. +* **`monoprop` has a mandatory Pauli-weight `cutoff`**, which the others have no + analogue of. The runner sets it to `n` — the whole register — so it never binds + and `lower_atol` is the only truncation in play. +* **PauliStrings.jl stores `im^{#Y}` inside the coefficient** (its `Matrix` + convention). The dump goes through `op_to_strings`, which puts its + coefficients on the same real footing as everyone else's. +* **PauliStrings.jl's own front door is `evolve(H, O, tspan; method=Trotter())`, + which we do not use** — it derives the gate list from the Hamiltonian's + internal string order, which is not the spec's order. The runner builds the + `TrotterGate` vector by hand instead. +* **Julia is timed after a warm-up run** so the reported time excludes JIT. +* **Circuit construction is outside the timed region for every engine.** + monoprop's `propagate` re-expands its gate list internally on each call, which + no flag hoists out; measured at 0.0–2.2 % of its total at these widths, so it + is left in rather than worked around. +* Everything is single-threaded (`julia -t1`; `pauli-prop` is single-threaded by + design; `monoprop` is capped as above; `ppvm` here uses no Rayon). + +### Where the two truncation rules actually diverge + +Worth writing down, because the intuitive guesses about this are wrong and the +scalar `observable` will mislead you about all of them. + +A Pauli rotation sends `Q → cos θ·Q + sin θ·(iPQ)`, so a word receives +contributions from **at most two** sources: itself, scaled by `cos`, and its +partner `PQ`, scaled by `sin`. `ppvm` thresholds that sum; `monoprop` thresholds +each contribution as it is emitted. Two consequences follow, and they set the +whole shape of the difference: + +* `monoprop` discards a child whenever `|c|·sin θ < atol` even though its parent + survives, so it loses every child of a parent in the band + `atol < |c| < atol/sin θ`. The band's width is `1/sin θ` — **set by the + rotation angle alone**, with no dependence on the shape of the coefficient + distribution. +* Because a sum of two terms can exceed its larger member by at most 2×, the + mass `ppvm` rescues is `O(atol)` per word. Both engines' errors are therefore + anchored to the same threshold and cannot be decoupled: drive `atol` low enough + that `ppvm` is accurate and `monoprop` is accurate too. + +So the gap needs three things *at once*: small `sin θ` to widen the band, depth +to populate it and compound the loss, and a scrambled instance so the loss shows +up in the answer instead of averaging out. Missing any one of them, the gap +collapses to 1.0–1.2× or reverses. Measured, at `n=8`: + +| circuit | `atol` | `monoprop` error / `ppvm` error | +|---|---|---| +| `scramble`, `dt=0.05` (`θ ≤ 0.1`), 3 200 gates | `1e-4` | **2.96×** (2.89–3.07× over 10 seeds) | +| `heisenberg`, `dt=0.05`, 320 steps | `1e-3` | 2.77× | +| `tfim`/`heisenberg`, `dt=0.1`, converged `atol` | `≤1e-5` | 1.02–1.15× | +| `scramble`, `dt=π/4` (random `θ ≤ π/2`), 80 gates | `1e-3…1e-5` | 0.85–1.00× | +| `heisenberg`, `dt=0.6` | `1e-3` | **0.67×** | + +[`accuracy_divergence.csv`](accuracy_divergence.csv) holds these sweeps, and +`plot_accuracy.py` renders them as four panels — the `atol` scaling, the depth +trend, the angle reversal, and the scrambler beside its control. + +The last two rows are the ones that catch people out. **Large angles favour +`monoprop`**, by up to 1.5×, because that is where `ppvm`'s rule has its own +failure mode: `truncate()` after each gate permanently deletes a term whose +merged coefficient has transiently cancelled below `atol`, where `monoprop` keeps +the row and lets later gates revive it. Uniform amplitudes do not help either — +`θ = π/4` is where the distribution is flattest (`sd(log|c|)` 0.66 against 2.23 +at `dt=0.1`) and it is also `monoprop`'s *best* regime, since `sin θ = 0.707` +makes the rejection band 1.41× wide, the narrowest possible. Flat amplitudes and +a wide band are mutually exclusive: you get the band from small angles, and small +angles force the geometric spread that makes amplitudes non-uniform. + +None of this is visible in the `observable` column. On the `scramble` instance +above the peak-error ratio has a median of 4.2× but a range of 0.1–576×, because +a single scalar's truncation error changes sign; the 576× is `ppvm` landing +accidentally near-exact on one seed. Neither engine loses the signal there — +relative errors are 0.4 % and 2.6 % against a peak of 0.16. + +## A run + +Apple M4 Pro (10 P + 4 E cores), macOS 26.6, rustc 1.96.0, Julia 1.12.6, +monoprop 0.8.0, PauliStrings.jl 1.10.1; single-threaded, `--steps 10 --dt 0.1 +--j 1 --h 1 --atol 1e-6 --iters 2`. Runtime relative to `ppvm`; **below 1.00× +means that engine beat `ppvm`**. `terms` is the shared support the three +accumulate-then-truncate engines all reach exactly. + +`tfim`, widths 8…64: + +| n | terms | ppvm | PauliPropagation.jl | PauliStrings.jl | pauli-prop | monoprop | +|---:|---:|---:|---:|---:|---:|---:| +| 8 | 4 701 | 0.0028 s | 2.42× | 4.90× | 11.13× | 0.98× | +| 16 | 19 529 | 0.0122 s | 3.10× | 6.26× | 5.53× | 0.69× | +| 24 | 34 353 | 0.0260 s | 5.78× | 8.05× | 5.32× | 0.54× | +| 32 | 49 177 | 0.0431 s | 4.98× | 8.39× | 7.05× | 0.45× | +| 40 | 64 001 | 0.0657 s | 6.33× | 9.20× | 5.28× | 0.43× | +| 48 | 78 825 | 0.0899 s | 7.32× | 9.02× | 5.37× | 0.40× | +| 56 | 93 649 | 0.1188 s | 10.46× | 9.84× | 5.33× | 0.35× | +| 64 | 108 473 | 0.1429 s | 11.26× | 13.18× | 5.71× | 0.33× | + +`heisenberg`, widths 6…14: + +| n | terms | ppvm | PauliPropagation.jl | PauliStrings.jl | pauli-prop | monoprop | +|---:|---:|---:|---:|---:|---:|---:| +| 6 | 1 022 | 0.0028 s | 4.28× | 6.77× | 5.02× | 1.03× | +| 8 | 16 324 | 0.0505 s | 3.43× | 6.26× | 2.02× | 0.83× | +| 10 | 225 353 | 0.5290 s | 3.11× | 6.79× | 1.51× | 0.77× | +| 12 | 1 174 849 | 2.7528 s | 2.59× | 6.50× | 1.13× | 0.67× | +| 14 | 2 915 879 | 6.0906 s | 2.23× | 7.09× | 0.88× | 0.54× | + +### Reading it + +**`ppvm` beats both Julia engines everywhere**, by 2.2–11.3× against +PauliPropagation.jl and 4.9–13.2× against PauliStrings.jl, and the margin widens +with `n` on TFIM. Those three carry an identical support term-for-term, so those +are clean ratios. + +**`ppvm` loses to `monoprop`** on both models — by up to 3.0× on TFIM at `n=64` +and 1.9× on Heisenberg at `n=14`, with the gap widening in `n` on both. On TFIM +the two carry the same support to within 1 %, so that column is a clean loss. On +Heisenberg `monoprop` is doing 4–12 % *more* above-threshold work than `ppvm` and +still finishing sooner, so the gap there is if anything understated — though it +is also tracking 2.3× as many rows in total, so the two engines are not making +the same space/time trade. + +**`pauli-prop` is 5.3–11.1× slower on TFIM but competitive on Heisenberg**, +reaching 0.88× at `n=14`. Read its column against the workload note: on TFIM at +`n=64` it is propagating 36 % fewer terms than the reference three, and on +Heisenberg `n=10`–`12` about 15 % more. + +## Scaling note + +The TFIM support grows roughly linearly in `n` at fixed depth, so that sweep +reaches 64 qubits cheaply. The Heisenberg support grows much faster — three +non-commuting bond rotations per bond per step — and the `Z₀` autocorrelator +itself converges once `n` exceeds the light cone (identical from `n≈10` at +`steps=10`), so widths past that measure scaling rather than new physics. + +## Files + +- `../../crates/ppvm-pauli-sum/examples/xbench.rs` — the `ppvm` runner. +- `../../julia-benchmarks/benches/xbench_pp.jl` — PauliPropagation.jl. +- `../../julia-benchmarks/benches/xbench_ps.jl` — PauliStrings.jl. +- `xbench_qiskit.py` — `pauli-prop`. +- `xbench_monoprop.py` — `monoprop`, including the thread-cap assertion. +- `run_xbench.py` — validation, driving, CSV merge, summary table. +- `plot_xbench.py` — the figure. +- `xbench_accuracy.py` — the sweeps behind `accuracy.csv` (all five engines vs + `atol`) and `accuracy_divergence.csv` (depth, angle, and the scrambler): each + engine's coefficient vector against a converged reference, so the timings can + be read next to what each truncation rule costs. Sweeps any of `--models`, + `--steps`, `--dt`, `--seeds`, `--atols` as comma-separated lists. +- `plot_accuracy.py` — the four-panel error-scaling figure from those two CSVs. + +[pp]: https://github.com/MSRudolph/PauliPropagation.jl +[ps]: https://github.com/nicolasloizeau/PauliStrings.jl +[qk]: https://github.com/Qiskit/pauli-prop +[mp]: https://github.com/Algorithmiq/monoprop diff --git a/benchmarks/cross-library/accuracy.csv b/benchmarks/cross-library/accuracy.csv new file mode 100644 index 000000000..4b2f3cbbd --- /dev/null +++ b/benchmarks/cross-library/accuracy.csv @@ -0,0 +1,51 @@ +model,library,qubits,steps,dt,seed,atol,ref_atol,ref_terms,ref_norm,terms,l2_err,l1_err,max_coeff_err,lost_mass,kept_subthreshold,dropped_above_atol,observable,ref_observable,obs_abs_err +tfim,ppvm,8,10,0.1,12345,1e-3,1e-16,27051,2.828427e+00,834,5.730178e-02,1.956886e+00,6.556158e-03,1.180394e+00,0,388,1.889668769617,1.888902654574,7.661150e-04 +tfim,pauli-propagation-jl,8,10,0.1,12345,1e-3,1e-16,27051,2.828427e+00,834,5.730178e-02,1.956886e+00,6.556158e-03,1.180394e+00,0,388,1.88966876961681,1.888902654574,7.661150e-04 +tfim,pauli-strings-jl,8,10,0.1,12345,1e-3,1e-16,27051,2.828427e+00,834,5.730178e-02,1.956886e+00,6.556158e-03,1.180394e+00,0,388,1.88966876961681,1.888902654574,7.661150e-04 +tfim,pauli-prop,8,10,0.1,12345,1e-3,1e-16,27051,2.828427e+00,743,9.747655e-02,3.201066e+00,8.889589e-03,1.519868e+00,0,479,1.9013036238629,1.888902654574,1.240097e-02 +tfim,monoprop,8,10,0.1,12345,1e-3,1e-16,27051,2.828427e+00,804,6.429517e-02,2.126366e+00,6.563651e-03,1.286381e+00,0,418,1.895370105228,1.888902654574,6.467451e-03 +tfim,ppvm,8,10,0.1,12345,1e-4,1e-16,27051,2.828427e+00,1730,7.713204e-03,3.399904e-01,5.538029e-04,1.657662e-01,0,558,1.887148304814,1.888902654574,1.754350e-03 +tfim,pauli-propagation-jl,8,10,0.1,12345,1e-4,1e-16,27051,2.828427e+00,1730,7.713204e-03,3.399904e-01,5.538029e-04,1.657662e-01,0,558,1.88714830481379,1.888902654574,1.754350e-03 +tfim,pauli-strings-jl,8,10,0.1,12345,1e-4,1e-16,27051,2.828427e+00,1730,7.713204e-03,3.399904e-01,5.538029e-04,1.657662e-01,0,558,1.88714830481379,1.888902654574,1.754350e-03 +tfim,pauli-prop,8,10,0.1,12345,1e-4,1e-16,27051,2.828427e+00,1527,1.655369e-02,6.909180e-01,1.396127e-03,2.516797e-01,0,761,1.89026934442961,1.888902654574,1.366690e-03 +tfim,monoprop,8,10,0.1,12345,1e-4,1e-16,27051,2.828427e+00,1705,8.349062e-03,3.553170e-01,6.351504e-04,1.787774e-01,1,584,1.888182135464,1.888902654574,7.205191e-04 +tfim,ppvm,8,10,0.1,12345,1e-5,1e-16,27051,2.828427e+00,3001,9.762681e-04,5.198656e-02,5.615238e-05,2.218986e-02,0,715,1.888799248154,1.888902654574,1.034064e-04 +tfim,pauli-propagation-jl,8,10,0.1,12345,1e-5,1e-16,27051,2.828427e+00,3001,9.762681e-04,5.198656e-02,5.615238e-05,2.218986e-02,0,715,1.88879924815445,1.888902654574,1.034064e-04 +tfim,pauli-strings-jl,8,10,0.1,12345,1e-5,1e-16,27051,2.828427e+00,3001,9.762681e-04,5.198656e-02,5.615238e-05,2.218986e-02,0,715,1.88879924815445,1.888902654574,1.034064e-04 +tfim,pauli-prop,8,10,0.1,12345,1e-5,1e-16,27051,2.828427e+00,2628,2.357978e-03,1.202946e-01,1.384898e-04,3.877258e-02,0,1088,1.88885571834255,1.888902654574,4.693623e-05 +tfim,monoprop,8,10,0.1,12345,1e-5,1e-16,27051,2.828427e+00,2998,1.023503e-03,5.342297e-02,5.628422e-05,2.245608e-02,0,718,1.888873511516,1.888902654574,2.914306e-05 +tfim,ppvm,8,10,0.1,12345,1e-6,1e-16,27051,2.828427e+00,4701,1.118273e-04,6.914840e-03,7.337651e-06,2.638321e-03,0,875,1.888902953057,1.888902654574,2.984830e-07 +tfim,pauli-propagation-jl,8,10,0.1,12345,1e-6,1e-16,27051,2.828427e+00,4701,1.118273e-04,6.914840e-03,7.337651e-06,2.638321e-03,0,875,1.88890295305684,1.888902654574,2.984828e-07 +tfim,pauli-strings-jl,8,10,0.1,12345,1e-6,1e-16,27051,2.828427e+00,4701,1.118273e-04,6.914840e-03,7.337651e-06,2.638321e-03,0,875,1.88890295305684,1.888902654574,2.984828e-07 +tfim,pauli-prop,8,10,0.1,12345,1e-6,1e-16,27051,2.828427e+00,4161,2.838847e-04,1.707832e-02,1.911237e-05,5.130690e-03,0,1415,1.88893188165528,1.888902654574,2.922708e-05 +tfim,monoprop,8,10,0.1,12345,1e-6,1e-16,27051,2.828427e+00,4662,1.182467e-04,7.215810e-03,7.389538e-06,2.782631e-03,1,915,1.888910840223,1.888902654574,8.185649e-06 +tfim,ppvm,8,10,0.1,12345,1e-7,1e-16,27051,2.828427e+00,6780,1.239548e-05,8.753233e-04,7.429950e-07,2.963253e-04,0,996,1.888901406583,1.888902654574,1.247991e-06 +tfim,pauli-propagation-jl,8,10,0.1,12345,1e-7,1e-16,27051,2.828427e+00,6780,1.239548e-05,8.753233e-04,7.429950e-07,2.963253e-04,0,996,1.88890140658276,1.888902654574,1.247991e-06 +tfim,pauli-strings-jl,8,10,0.1,12345,1e-7,1e-16,27051,2.828427e+00,6780,1.239548e-05,8.753233e-04,7.429950e-07,2.963253e-04,0,996,1.88890140658276,1.888902654574,1.247991e-06 +tfim,pauli-prop,8,10,0.1,12345,1e-7,1e-16,27051,2.828427e+00,6056,3.227852e-05,2.215272e-03,1.725103e-06,6.248028e-04,0,1720,1.88890165053338,1.888902654574,1.004041e-06 +tfim,monoprop,8,10,0.1,12345,1e-7,1e-16,27051,2.828427e+00,6753,1.309534e-05,9.190766e-04,7.516440e-07,3.063929e-04,1,1024,1.888902140509,1.888902654574,5.140650e-07 +heisenberg,ppvm,8,10,0.1,12345,1e-3,1e-16,16384,1.000000e+00,4566,1.935867e-01,1.683728e+01,6.612049e-03,8.467940e+00,370,2994,0.1843817263679,0.1835878942086,7.938322e-04 +heisenberg,pauli-propagation-jl,8,10,0.1,12345,1e-3,1e-16,16384,1.000000e+00,4566,1.935867e-01,1.683728e+01,6.612049e-03,8.467940e+00,370,2994,0.184381726367865,0.1835878942086,7.938322e-04 +heisenberg,pauli-strings-jl,8,10,0.1,12345,1e-3,1e-16,16384,1.000000e+00,4566,1.935867e-01,1.683728e+01,6.612049e-03,8.467940e+00,370,2994,0.184381726367865,0.1835878942086,7.938322e-04 +heisenberg,pauli-prop,8,10,0.1,12345,1e-3,1e-16,16384,1.000000e+00,5072,2.384650e-01,2.082594e+01,1.022082e-02,9.244654e+00,1051,3169,0.186064764275212,0.1835878942086,2.476870e-03 +heisenberg,monoprop,8,10,0.1,12345,1e-3,1e-16,16384,1.000000e+00,5106,1.937250e-01,1.714593e+01,6.293714e-03,8.315984e+00,836,2920,0.18421262641,0.1835878942086,6.247322e-04 +heisenberg,ppvm,8,10,0.1,12345,1e-4,1e-16,16384,1.000000e+00,11923,2.564185e-02,2.511013e+00,1.063198e-03,5.963471e-01,1107,2274,0.1834786826111,0.1835878942086,1.092116e-04 +heisenberg,pauli-propagation-jl,8,10,0.1,12345,1e-4,1e-16,16384,1.000000e+00,11923,2.564185e-02,2.511013e+00,1.063198e-03,5.963471e-01,1107,2274,0.183478682611051,0.1835878942086,1.092116e-04 +heisenberg,pauli-strings-jl,8,10,0.1,12345,1e-4,1e-16,16384,1.000000e+00,11923,2.564185e-02,2.511013e+00,1.063198e-03,5.963471e-01,1107,2274,0.183478682611051,0.1835878942086,1.092116e-04 +heisenberg,pauli-prop,8,10,0.1,12345,1e-4,1e-16,16384,1.000000e+00,13628,4.451207e-02,4.398424e+00,2.046005e-03,6.386030e-01,2836,2298,0.183655881056651,0.1835878942086,6.798685e-05 +heisenberg,monoprop,8,10,0.1,12345,1e-4,1e-16,16384,1.000000e+00,13143,2.725847e-02,2.714465e+00,1.193874e-03,5.551494e-01,2074,2021,0.183546769437,0.1835878942086,4.112477e-05 +heisenberg,ppvm,8,10,0.1,12345,1e-5,1e-16,16384,1.000000e+00,16101,2.426734e-03,2.449232e-01,1.074420e-04,1.125860e-02,963,420,0.1835902065899,0.1835878942086,2.312381e-06 +heisenberg,pauli-propagation-jl,8,10,0.1,12345,1e-5,1e-16,16384,1.000000e+00,16101,2.426734e-03,2.449232e-01,1.074420e-04,1.125860e-02,963,420,0.18359020658986,0.1835878942086,2.312381e-06 +heisenberg,pauli-strings-jl,8,10,0.1,12345,1e-5,1e-16,16384,1.000000e+00,16101,2.426734e-03,2.449232e-01,1.074420e-04,1.125860e-02,963,420,0.183590206589861,0.1835878942086,2.312381e-06 +heisenberg,pauli-prop,8,10,0.1,12345,1e-5,1e-16,16384,1.000000e+00,17803,4.507725e-03,4.730143e-01,1.675695e-04,1.250100e-02,2677,432,0.183584912095183,0.1835878942086,2.982113e-06 +heisenberg,monoprop,8,10,0.1,12345,1e-5,1e-16,16384,1.000000e+00,17223,2.557524e-03,2.663380e-01,1.096506e-04,9.925007e-03,2027,362,0.183585575353,0.1835878942086,2.318856e-06 +heisenberg,ppvm,8,10,0.1,12345,1e-6,1e-16,16384,1.000000e+00,16324,1.873617e-04,1.844112e-02,8.671387e-06,1.957258e-04,191,82,0.1835880798476,0.1835878942086,1.856390e-07 +heisenberg,pauli-propagation-jl,8,10,0.1,12345,1e-6,1e-16,16384,1.000000e+00,16324,1.873617e-04,1.844112e-02,8.671387e-06,1.957258e-04,191,82,0.183588079847589,0.1835878942086,1.856390e-07 +heisenberg,pauli-strings-jl,8,10,0.1,12345,1e-6,1e-16,16384,1.000000e+00,16324,1.873617e-04,1.844112e-02,8.671387e-06,1.957258e-04,191,82,0.183588079847589,0.1835878942086,1.856390e-07 +heisenberg,pauli-prop,8,10,0.1,12345,1e-6,1e-16,16384,1.000000e+00,17004,3.764296e-04,3.802784e-02,1.821500e-05,2.767212e-04,894,105,0.183587851042604,0.1835878942086,4.316600e-08 +heisenberg,monoprop,8,10,0.1,12345,1e-6,1e-16,16384,1.000000e+00,16910,1.993482e-04,2.007446e-02,1.047708e-05,1.937806e-04,774,79,0.183587892322,0.1835878942086,1.886600e-09 +heisenberg,ppvm,8,10,0.1,12345,1e-7,1e-16,16384,1.000000e+00,16377,1.541357e-05,1.514858e-03,7.285205e-07,3.268561e-06,39,14,0.1835878872944,0.1835878942086,6.914200e-09 +heisenberg,pauli-propagation-jl,8,10,0.1,12345,1e-7,1e-16,16384,1.000000e+00,16377,1.541357e-05,1.514858e-03,7.285205e-07,3.268561e-06,39,14,0.183587887294405,0.1835878942086,6.914195e-09 +heisenberg,pauli-strings-jl,8,10,0.1,12345,1e-7,1e-16,16384,1.000000e+00,16377,1.541357e-05,1.514858e-03,7.285205e-07,3.268561e-06,39,14,0.183587887294405,0.1835878942086,6.914195e-09 +heisenberg,pauli-prop,8,10,0.1,12345,1e-7,1e-16,16384,1.000000e+00,16528,3.026389e-05,3.008580e-03,1.430835e-06,5.571167e-06,195,19,0.18358789229432,0.1835878942086,1.914280e-09 +heisenberg,monoprop,8,10,0.1,12345,1e-7,1e-16,16384,1.000000e+00,16558,1.539943e-05,1.536256e-03,8.884821e-07,3.199623e-06,217,11,0.183587893774,0.1835878942086,4.346000e-10 diff --git a/benchmarks/cross-library/accuracy_divergence.csv b/benchmarks/cross-library/accuracy_divergence.csv new file mode 100644 index 000000000..669426026 --- /dev/null +++ b/benchmarks/cross-library/accuracy_divergence.csv @@ -0,0 +1,99 @@ +model,library,qubits,steps,dt,seed,atol,ref_atol,ref_terms,ref_norm,terms,l2_err,l1_err,max_coeff_err,lost_mass,kept_subthreshold,dropped_above_atol,observable,ref_observable,obs_abs_err +heisenberg,ppvm,6,10,0.05,12345,1e-3,1e-16,1024,1.000000e+00,222,4.036339e-02,8.303254e-01,3.634578e-03,4.669530e-01,0,155,0.3832181169286,0.3822563306389,9.617863e-04 +heisenberg,monoprop,6,10,0.05,12345,1e-3,1e-16,1024,1.000000e+00,224,4.214458e-02,8.422781e-01,4.098614e-03,4.650453e-01,2,155,0.383158101684,0.3822563306389,9.017710e-04 +heisenberg,ppvm,6,10,0.05,12345,1e-5,1e-16,1024,1.000000e+00,826,5.404658e-04,1.318947e-02,6.420822e-05,1.748059e-03,18,76,0.3822591781592,0.3822563306389,2.847520e-06 +heisenberg,monoprop,6,10,0.05,12345,1e-5,1e-16,1024,1.000000e+00,906,5.160546e-04,1.273802e-02,6.591156e-05,1.255044e-03,74,52,0.382257798983,0.3822563306389,1.468344e-06 +heisenberg,ppvm,6,20,0.05,12345,1e-3,1e-16,1024,1.000000e+00,866,1.070213e-01,2.716866e+00,1.192541e-02,4.248833e-01,90,130,0.1859379068078,0.1853914774626,5.464293e-04 +heisenberg,monoprop,6,20,0.05,12345,1e-3,1e-16,1024,1.000000e+00,993,1.252160e-01,3.207887e+00,1.363090e-02,3.507914e-01,180,93,0.185871276693,0.1853914774626,4.797992e-04 +heisenberg,ppvm,6,20,0.05,12345,1e-5,1e-16,1024,1.000000e+00,1039,6.429122e-04,1.651200e-02,7.105414e-05,1.100133e-04,20,3,0.1854168229734,0.1853914774626,2.534551e-05 +heisenberg,monoprop,6,20,0.05,12345,1e-5,1e-16,1024,1.000000e+00,1167,7.136553e-04,1.958904e-02,7.192523e-05,1.181164e-04,149,4,0.185416479607,0.1853914774626,2.500214e-05 +heisenberg,ppvm,6,40,0.05,12345,1e-3,1e-16,1024,1.000000e+00,1177,1.474150e-01,4.066484e+00,1.613299e-02,2.478436e-01,239,49,0.1284102386314,0.1308659700642,2.455731e-03 +heisenberg,monoprop,6,40,0.05,12345,1e-3,1e-16,1024,1.000000e+00,1430,2.353190e-01,7.134715e+00,2.271986e-02,2.480553e-01,486,43,0.130211985822,0.1308659700642,6.539842e-04 +heisenberg,ppvm,6,40,0.05,12345,1e-5,1e-16,1024,1.000000e+00,1024,6.308929e-04,1.575149e-02,7.857540e-05,0,0,0,0.130870034274,0.1308659700642,4.064210e-06 +heisenberg,monoprop,6,40,0.05,12345,1e-5,1e-16,1024,1.000000e+00,1099,6.908877e-04,1.833722e-02,6.336108e-05,0,75,0,0.130868964327,0.1308659700642,2.994263e-06 +heisenberg,ppvm,6,80,0.05,12345,1e-3,1e-16,1024,1.000000e+00,1236,1.646609e-01,4.623798e+00,2.112876e-02,1.582509e-01,268,34,0.08887939929655,0.09030014507815,1.420746e-03 +heisenberg,monoprop,6,80,0.05,12345,1e-3,1e-16,1024,1.000000e+00,1638,3.087508e-01,9.415847e+00,3.496292e-02,2.588310e-01,675,39,0.075129187818,0.09030014507815,1.517096e-02 +heisenberg,ppvm,6,80,0.05,12345,1e-5,1e-16,1024,1.000000e+00,1024,6.290098e-04,1.561913e-02,7.908085e-05,0,0,0,0.09032734267639,0.09030014507815,2.719760e-05 +heisenberg,monoprop,6,80,0.05,12345,1e-5,1e-16,1024,1.000000e+00,1035,6.629870e-04,1.658641e-02,7.977269e-05,0,11,0,0.090328492789,0.09030014507815,2.834771e-05 +heisenberg,ppvm,6,160,0.05,12345,1e-3,1e-16,1024,1.000000e+00,1213,1.813172e-01,5.016238e+00,2.172605e-02,1.693484e-01,246,25,0.09555696196697,0.09803624235248,2.479280e-03 +heisenberg,monoprop,6,160,0.05,12345,1e-3,1e-16,1024,1.000000e+00,1691,4.000477e-01,1.195438e+01,4.531704e-02,3.583095e-01,735,36,0.118446290328,0.09803624235248,2.041005e-02 +heisenberg,ppvm,6,160,0.05,12345,1e-5,1e-16,1024,1.000000e+00,1024,6.334569e-04,1.596263e-02,6.675496e-05,0,0,0,0.09802627447504,0.09803624235248,9.967877e-06 +heisenberg,monoprop,6,160,0.05,12345,1e-5,1e-16,1024,1.000000e+00,1022,6.520457e-04,1.635314e-02,7.891152e-05,7.021406e-05,0,2,0.098031933021,0.09803624235248,4.309331e-06 +heisenberg,ppvm,6,320,0.05,12345,1e-3,1e-16,1024,1.000000e+00,1235,2.375216e-01,6.534776e+00,2.462977e-02,3.085507e-01,275,46,0.2136063818482,0.2171604004367,3.554019e-03 +heisenberg,monoprop,6,320,0.05,12345,1e-3,1e-16,1024,1.000000e+00,1787,6.576496e-01,1.928416e+01,7.418257e-02,7.182283e-01,821,40,0.192691392802,0.2171604004367,2.446901e-02 +heisenberg,ppvm,6,320,0.05,12345,1e-5,1e-16,1024,1.000000e+00,1024,6.369157e-04,1.611839e-02,7.839926e-05,0,0,0,0.2171604704937,0.2171604004367,7.005700e-08 +heisenberg,monoprop,6,320,0.05,12345,1e-5,1e-16,1024,1.000000e+00,1024,6.487131e-04,1.640341e-02,7.181142e-05,0,0,0,0.217165017466,0.2171604004367,4.617029e-06 +heisenberg,ppvm,8,10,0.02,12345,1e-3,1e-16,16561,1.000000e+00,27,1.365602e-02,1.124382e-01,4.679180e-03,6.805587e-02,0,26,0.8522590829718,0.8522087425321,5.034044e-05 +heisenberg,monoprop,8,10,0.02,12345,1e-3,1e-16,16561,1.000000e+00,25,1.452084e-02,1.197392e-01,5.736855e-03,7.172758e-02,0,28,0.852274606892,0.8522087425321,6.586436e-05 +heisenberg,ppvm,8,10,0.05,12345,1e-3,1e-16,16384,1.000000e+00,222,4.577082e-02,1.518473e+00,3.673186e-03,1.197095e+00,0,263,0.3832181169286,0.3822563787026,9.617382e-04 +heisenberg,monoprop,8,10,0.05,12345,1e-3,1e-16,16384,1.000000e+00,224,4.684098e-02,1.525291e+00,3.921037e-03,1.195189e+00,2,263,0.383158101684,0.3822563787026,9.017230e-04 +heisenberg,ppvm,8,10,0.1,12345,1e-3,1e-16,16384,1.000000e+00,4566,1.935867e-01,1.683728e+01,6.612049e-03,8.467940e+00,370,2994,0.1843817263679,0.1835878942086,7.938322e-04 +heisenberg,monoprop,8,10,0.1,12345,1e-3,1e-16,16384,1.000000e+00,5106,1.937250e-01,1.714593e+01,6.293714e-03,8.315984e+00,836,2920,0.18421262641,0.1835878942086,6.247322e-04 +heisenberg,ppvm,8,10,0.2,12345,1e-3,1e-16,16384,1.000000e+00,9617,3.392238e-01,3.496440e+01,1.721673e-02,1.503669e+01,2670,5303,0.1206720275243,0.1197752470608,8.967805e-04 +heisenberg,monoprop,8,10,0.2,12345,1e-3,1e-16,16384,1.000000e+00,9516,3.656091e-01,3.825462e+01,1.981698e-02,1.708788e+01,3084,5818,0.119257133464,0.1197752470608,5.181136e-04 +heisenberg,ppvm,8,10,0.4,12345,1e-3,1e-16,16384,1.000000e+00,12275,3.954437e-01,4.123504e+01,1.592698e-02,1.312605e+01,2316,3979,0.06411724002328,0.06498038532114,8.631453e-04 +heisenberg,monoprop,8,10,0.4,12345,1e-3,1e-16,16384,1.000000e+00,11524,3.180263e-01,3.235889e+01,1.186265e-02,1.048401e+01,1169,3583,0.064849564292,0.06498038532114,1.308210e-04 +heisenberg,ppvm,8,10,0.6,12345,1e-3,1e-16,16384,1.000000e+00,13027,3.794862e-01,3.972636e+01,1.750531e-02,1.173227e+01,3457,3830,0.1472159226886,0.1492943370647,2.078414e-03 +heisenberg,monoprop,8,10,0.6,12345,1e-3,1e-16,16384,1.000000e+00,12431,2.526846e-01,2.554907e+01,1.121828e-02,6.632730e+00,1625,2594,0.149259337365,0.1492943370647,3.499970e-05 +heisenberg,ppvm,8,10,0.785,12345,1e-3,1e-16,7322,1.000000e+00,1,9.422127e-03,3.564081e-02,3.185207e-03,3.561038e-02,0,14,0,5.707176610115e-06,5.707177e-06 +heisenberg,monoprop,8,10,0.785,12345,1e-3,1e-16,7322,1.000000e+00,1,9.422127e-03,3.564081e-02,3.185207e-03,3.561038e-02,0,14,0,5.707176610115e-06,5.707177e-06 +scramble,ppvm,8,400,0.05,1,1e-3,1e-16,65535,1.000000e+00,7145,9.460647e-01,1.886151e+02,2.247844e-02,1.588307e+02,1273,44679,0.1790007556962,0.1614730211106,1.752773e-02 +scramble,monoprop,8,400,0.05,1,1e-3,1e-16,65535,1.000000e+00,7923,9.619359e-01,1.897827e+02,2.840829e-02,1.554452e+02,1401,44029,0.167740351377,0.1614730211106,6.267330e-03 +scramble,ppvm,8,400,0.05,1,1e-4,1e-16,65535,1.000000e+00,64007,1.366369e-01,2.789267e+01,2.364321e-03,6.780471e-01,1332,1293,0.1614802576756,0.1614730211106,7.236565e-06 +scramble,monoprop,8,400,0.05,1,1e-4,1e-16,65535,1.000000e+00,63920,4.179348e-01,8.505128e+01,9.345567e-03,2.067460e+00,1480,1528,0.165644417085,0.1614730211106,4.171396e-03 +scramble,ppvm,8,400,0.05,1,1e-5,1e-16,65535,1.000000e+00,65388,7.801179e-03,1.593150e+00,1.427393e-04,3.874730e-03,110,119,0.1614469519445,0.1614730211106,2.606917e-05 +scramble,monoprop,8,400,0.05,1,1e-5,1e-16,65535,1.000000e+00,65376,1.940221e-02,3.961147e+00,3.854044e-04,9.243616e-03,122,143,0.161298097215,0.1614730211106,1.749239e-04 +scramble,ppvm,8,400,0.05,2,1e-3,1e-16,65535,1.000000e+00,7277,9.563637e-01,1.915102e+02,2.113910e-02,1.610994e+02,1276,44811,0.145258196215,0.141763190312,3.495006e-03 +scramble,monoprop,8,400,0.05,2,1e-3,1e-16,65535,1.000000e+00,7884,9.710104e-01,1.920922e+02,2.926046e-02,1.580419e+02,1338,44266,0.121961388302,0.141763190312,1.980180e-02 +scramble,ppvm,8,400,0.05,2,1e-4,1e-16,65535,1.000000e+00,64004,1.406716e-01,2.866367e+01,2.309150e-03,6.676462e-01,1275,1297,0.1411270619147,0.141763190312,6.361284e-04 +scramble,monoprop,8,400,0.05,2,1e-4,1e-16,65535,1.000000e+00,63950,4.073952e-01,8.292705e+01,7.614204e-03,1.971686e+00,1420,1496,0.138618759706,0.141763190312,3.144431e-03 +scramble,ppvm,8,400,0.05,2,1e-5,1e-16,65535,1.000000e+00,65400,8.667319e-03,1.771308e+00,1.364485e-04,3.740238e-03,106,96,0.1417936293827,0.141763190312,3.043907e-05 +scramble,monoprop,8,400,0.05,2,1e-5,1e-16,65535,1.000000e+00,65398,1.939419e-02,3.962640e+00,3.397129e-04,8.654144e-03,124,116,0.141898202448,0.141763190312,1.350121e-04 +scramble,ppvm,8,400,0.05,3,1e-3,1e-16,65535,1.000000e+00,7384,9.450019e-01,1.887800e+02,2.825414e-02,1.580012e+02,1296,44432,0.1524323856682,0.1460185792372,6.413806e-03 +scramble,monoprop,8,400,0.05,3,1e-3,1e-16,65535,1.000000e+00,8026,9.643916e-01,1.901604e+02,3.122063e-02,1.550255e+02,1377,43871,0.142537145171,0.1460185792372,3.481434e-03 +scramble,ppvm,8,400,0.05,3,1e-4,1e-16,65535,1.000000e+00,63904,1.404446e-01,2.871432e+01,2.399302e-03,7.247692e-01,1277,1410,0.1452846799804,0.1460185792372,7.338993e-04 +scramble,monoprop,8,400,0.05,3,1e-4,1e-16,65535,1.000000e+00,64008,4.156444e-01,8.458307e+01,8.494982e-03,1.963916e+00,1432,1461,0.140959989003,0.1460185792372,5.058590e-03 +scramble,ppvm,8,400,0.05,3,1e-5,1e-16,65535,1.000000e+00,65376,8.038709e-03,1.642135e+00,1.416798e-04,4.002692e-03,101,117,0.1460074411703,0.1460185792372,1.113807e-05 +scramble,monoprop,8,400,0.05,3,1e-5,1e-16,65535,1.000000e+00,65378,1.951992e-02,3.980961e+00,3.259493e-04,9.704331e-03,124,138,0.146049390259,0.1460185792372,3.081102e-05 +scramble,ppvm,8,400,0.05,4,1e-3,1e-16,65535,1.000000e+00,7409,9.518168e-01,1.894747e+02,2.758227e-02,1.578000e+02,1251,44408,0.1682610320682,0.1693302302423,1.069198e-03 +scramble,monoprop,8,400,0.05,4,1e-3,1e-16,65535,1.000000e+00,8156,9.682760e-01,1.903117e+02,4.164777e-02,1.539202e+02,1307,43717,0.144043912834,0.1693302302423,2.528632e-02 +scramble,ppvm,8,400,0.05,4,1e-4,1e-16,65535,1.000000e+00,63869,1.409938e-01,2.882227e+01,2.557839e-03,7.418089e-01,1298,1408,0.1686796776657,0.1693302302423,6.505526e-04 +scramble,monoprop,8,400,0.05,4,1e-4,1e-16,65535,1.000000e+00,63968,4.106727e-01,8.365478e+01,7.280307e-03,1.928183e+00,1464,1475,0.165243860705,0.1693302302423,4.086370e-03 +scramble,ppvm,8,400,0.05,4,1e-5,1e-16,65535,1.000000e+00,65376,8.336955e-03,1.701727e+00,1.487382e-04,4.134139e-03,123,126,0.1693221010763,0.1693302302423,8.129166e-06 +scramble,monoprop,8,400,0.05,4,1e-5,1e-16,65535,1.000000e+00,65368,2.026217e-02,4.130641e+00,3.527960e-04,1.150466e-02,141,152,0.169290141837,0.1693302302423,4.008841e-05 +scramble,ppvm,8,400,0.05,5,1e-3,1e-16,65535,1.000000e+00,7168,9.504922e-01,1.888251e+02,2.381869e-02,1.581427e+02,1314,44541,0.1488295179599,0.1622150514644,1.338553e-02 +scramble,monoprop,8,400,0.05,5,1e-3,1e-16,65535,1.000000e+00,8263,9.733226e-01,1.903354e+02,4.313585e-02,1.532102e+02,1398,43530,0.187629107769,0.1622150514644,2.541406e-02 +scramble,ppvm,8,400,0.05,5,1e-4,1e-16,65535,1.000000e+00,63819,1.413663e-01,2.886298e+01,2.541082e-03,7.565301e-01,1294,1459,0.1601476023761,0.1622150514644,2.067449e-03 +scramble,monoprop,8,400,0.05,5,1e-4,1e-16,65535,1.000000e+00,63888,4.233586e-01,8.599954e+01,8.111627e-03,2.104861e+00,1458,1554,0.154900547436,0.1622150514644,7.314504e-03 +scramble,ppvm,8,400,0.05,5,1e-5,1e-16,65535,1.000000e+00,65372,7.935054e-03,1.618603e+00,1.321421e-04,4.264874e-03,115,128,0.1621842282147,0.1622150514644,3.082325e-05 +scramble,monoprop,8,400,0.05,5,1e-5,1e-16,65535,1.000000e+00,65383,1.933189e-02,3.942074e+00,3.162846e-04,9.839686e-03,141,143,0.162116789286,0.1622150514644,9.826218e-05 +scramble,ppvm,8,10,0.7853981633974483,1,1e-3,1e-16,65535,1.000000e+00,50657,3.090565e-01,6.278704e+01,5.755294e-03,1.610769e+01,5920,6774,0.002947642913682,0.004567937262123,1.620294e-03 +scramble,monoprop,8,10,0.7853981633974483,1,1e-3,1e-16,65535,1.000000e+00,49970,2.637431e-01,5.349200e+01,5.279884e-03,1.500207e+01,4766,6307,0.002980535216,0.004567937262123,1.587402e-03 +scramble,ppvm,8,10,0.7853981633974483,1,1e-4,1e-16,65535,1.000000e+00,64150,1.674435e-02,3.406150e+00,3.050943e-04,9.712100e-02,374,342,0.004554433189204,0.004567937262123,1.350407e-05 +scramble,monoprop,8,10,0.7853981633974483,1,1e-4,1e-16,65535,1.000000e+00,64144,1.270353e-02,2.583770e+00,2.536806e-04,8.731900e-02,294,268,0.004512298994,0.004567937262123,5.563827e-05 +scramble,ppvm,8,10,0.7853981633974483,1,1e-5,1e-16,65535,1.000000e+00,65396,1.041360e-03,2.120041e-01,2.233626e-05,7.722638e-04,26,21,0.004568585960411,0.004567937262123,6.486983e-07 +scramble,monoprop,8,10,0.7853981633974483,1,1e-5,1e-16,65535,1.000000e+00,65392,9.593503e-04,1.913561e-01,2.030360e-05,7.639987e-04,18,17,0.004566272673,0.004567937262123,1.664589e-06 +scramble,ppvm,8,10,0.7853981633974483,2,1e-3,1e-16,65535,1.000000e+00,48647,2.791048e-01,5.628207e+01,5.037719e-03,1.722668e+01,5167,7348,-0.003044049375097,-0.002123511847918,9.205375e-04 +scramble,monoprop,8,10,0.7853981633974483,2,1e-3,1e-16,65535,1.000000e+00,49443,2.734635e-01,5.507724e+01,5.526536e-03,1.610547e+01,5454,6839,-0.002009768618,-0.002123511847918,1.137432e-04 +scramble,ppvm,8,10,0.7853981633974483,2,1e-4,1e-16,65535,1.000000e+00,64060,1.508072e-02,3.042391e+00,3.327267e-04,9.911472e-02,394,348,-0.002128546499306,-0.002123511847918,5.034651e-06 +scramble,monoprop,8,10,0.7853981633974483,2,1e-4,1e-16,65535,1.000000e+00,64063,1.363255e-02,2.758181e+00,3.033626e-04,9.342819e-02,353,304,-0.002116198144,-0.002123511847918,7.313704e-06 +scramble,ppvm,8,10,0.7853981633974483,2,1e-5,1e-16,65535,1.000000e+00,65377,8.178769e-04,1.652772e-01,1.911102e-05,8.209238e-04,27,15,-0.002124049692306,-0.002123511847918,5.378444e-07 +scramble,monoprop,8,10,0.7853981633974483,2,1e-5,1e-16,65535,1.000000e+00,65362,6.915502e-04,1.389042e-01,2.041713e-05,9.752001e-04,18,21,-0.002124065586,-0.002123511847918,5.537381e-07 +scramble,ppvm,8,10,0.7853981633974483,3,1e-3,1e-16,65534,1.000000e+00,23800,1.482153e-01,2.838560e+01,3.781926e-03,1.934654e+01,2061,4604,0.01021329788007,0.01011042655569,1.028713e-04 +scramble,monoprop,8,10,0.7853981633974483,3,1e-3,1e-16,65534,1.000000e+00,23958,1.416023e-01,2.720370e+01,3.105871e-03,1.887551e+01,1907,4292,0.010168429392,0.01011042655569,5.800284e-05 +scramble,ppvm,8,10,0.7853981633974483,3,1e-4,1e-16,65534,1.000000e+00,57545,1.677474e-02,3.310245e+00,4.291547e-04,5.726769e-01,1649,2112,0.01012727159909,0.01011042655569,1.684504e-05 +scramble,monoprop,8,10,0.7853981633974483,3,1e-4,1e-16,65534,1.000000e+00,57569,1.497634e-02,2.933245e+00,4.037880e-04,5.407710e-01,1440,1879,0.010170823571,0.01011042655569,6.039702e-05 +scramble,ppvm,8,10,0.7853981633974483,3,1e-5,1e-16,65534,1.000000e+00,64720,1.185833e-03,2.343507e-01,2.743055e-05,5.164075e-03,175,166,0.01011032002169,0.01011042655569,1.065340e-07 +scramble,monoprop,8,10,0.7853981633974483,3,1e-5,1e-16,65534,1.000000e+00,64695,9.601826e-04,1.883416e-01,2.589419e-05,5.002271e-03,133,149,0.010108359939,0.01011042655569,2.066617e-06 +scramble,ppvm,8,10,0.7853981633974483,4,1e-3,1e-16,65535,1.000000e+00,44916,2.344832e-01,4.721697e+01,4.659655e-03,1.772357e+01,4861,7238,-0.002857638598484,-0.002233897581292,6.237410e-04 +scramble,monoprop,8,10,0.7853981633974483,4,1e-3,1e-16,65535,1.000000e+00,44165,2.366092e-01,4.750515e+01,4.174483e-03,1.888975e+01,4683,7811,-0.002729812368,-0.002233897581292,4.959148e-04 +scramble,ppvm,8,10,0.7853981633974483,4,1e-4,1e-16,65535,1.000000e+00,63664,1.656759e-02,3.339823e+00,3.475736e-04,1.364974e-01,544,528,-0.002347368442479,-0.002233897581292,1.134709e-04 +scramble,monoprop,8,10,0.7853981633974483,4,1e-4,1e-16,65535,1.000000e+00,63678,1.610363e-02,3.253601e+00,3.508339e-04,1.298318e-01,494,464,-0.002305939048,-0.002233897581292,7.204147e-05 +scramble,ppvm,8,10,0.7853981633974483,4,1e-5,1e-16,65535,1.000000e+00,65344,1.076490e-03,2.173143e-01,2.137590e-05,1.167930e-03,44,36,-0.0022253215307,-0.002233897581292,8.576051e-06 +scramble,monoprop,8,10,0.7853981633974483,4,1e-5,1e-16,65535,1.000000e+00,65344,9.685830e-04,1.950917e-01,2.328722e-05,1.113917e-03,37,29,-0.002225810421,-0.002233897581292,8.087160e-06 +scramble,ppvm,8,10,0.7853981633974483,5,1e-3,1e-16,65535,1.000000e+00,46519,2.524708e-01,5.087362e+01,4.986171e-03,1.754877e+01,5039,7322,0,-0.001223130468711,1.223130e-03 +scramble,monoprop,8,10,0.7853981633974483,5,1e-3,1e-16,65535,1.000000e+00,47075,2.425934e-01,4.885867e+01,4.802603e-03,1.636502e+01,5028,6755,0,-0.001223130468711,1.223130e-03 +scramble,ppvm,8,10,0.7853981633974483,5,1e-4,1e-16,65535,1.000000e+00,63800,1.575686e-02,3.174853e+00,3.102006e-04,1.252977e-01,441,469,-0.001210275579338,-0.001223130468711,1.285489e-05 +scramble,monoprop,8,10,0.7853981633974483,5,1e-4,1e-16,65535,1.000000e+00,63817,1.319639e-02,2.629690e+00,2.919619e-04,1.111590e-01,347,358,-0.001235391553,-0.001223130468711,1.226108e-05 +scramble,ppvm,8,10,0.7853981633974483,5,1e-5,1e-16,65535,1.000000e+00,65366,8.791609e-04,1.777130e-01,2.111646e-05,8.729767e-04,28,16,-0.001221170841282,-0.001223130468711,1.959627e-06 +scramble,monoprop,8,10,0.7853981633974483,5,1e-5,1e-16,65535,1.000000e+00,65365,6.298597e-04,1.269766e-01,1.810326e-05,8.522953e-04,22,11,-0.001224291877,-0.001223130468711,1.161408e-06 diff --git a/benchmarks/cross-library/plot_accuracy.py b/benchmarks/cross-library/plot_accuracy.py new file mode 100644 index 000000000..0ec9135ec --- /dev/null +++ b/benchmarks/cross-library/plot_accuracy.py @@ -0,0 +1,315 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 +"""Render the truncation-error scaling from `xbench_accuracy.py`'s CSVs. + +Four panels, because four separate questions were asked of this data and each one +has a different answer: + +1. **Error vs `atol`** — every engine falls one decade per decade of `atol`, so + the two truncation rules share a convergence order and differ only in the + constant. This panel is why the rules are not qualitatively apart. +2. **Error ratio vs depth** — the loss compounds per gate, so the gap grows with + circuit depth. +3. **Error ratio vs rotation angle** — the ratio *crosses 1.0*. Large angles + favour monoprop, because that is where accumulate-then-truncate has its own + failure mode: `truncate()` after each gate permanently deletes a term whose + merged coefficient transiently cancelled. Any summary quoting one ratio is + hiding this panel. +4. **The scrambling instance** — where the rules genuinely separate, beside a + `θ = π/4` control on the same circuit family that reverses the sign. + +Ratios are monoprop / ppvm on the `L2` distance to a converged reference, so +above 1.0 means monoprop is less accurate. The scalar `observable` is +deliberately not plotted: its truncation error changes sign, so ratios built from +it are dominated by whichever engine happened to land near a zero crossing (see +`README.md`). + + uv run --no-project --with matplotlib python3 plot_accuracy.py \ + --csv accuracy.csv --divergence accuracy_divergence.csv \ + --out target/xbench/accuracy.png +""" + +from __future__ import annotations + +import argparse +import csv +import statistics +import textwrap +from collections import defaultdict +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from plot_xbench import SERIES, SURFACE, TEXT_PRIMARY, TEXT_SECONDARY, style + +BASELINE = "ppvm" +PAIR = ("ppvm", "monoprop") +# Ratio panels plot a *derived* quantity, not an engine, so they stay off the +# engine palette — those hues are reserved for identity throughout the harness. +RATIO_INK = ["#0b0b0b", "#8a8986"] +RATIO_MARKS = ["v", "o"] + + +def load(path: Path) -> list[dict[str, str]]: + with path.open() as fh: + return list(csv.DictReader(fh)) + + +def relative(rows, library: str) -> dict[float, list[float]]: + """`{atol: [relative L2 error, ...]}` for one engine, pooled over seeds.""" + out: dict[float, list[float]] = defaultdict(list) + for r in rows: + if r["library"] == library: + out[float(r["atol"])].append(float(r["l2_err"]) / float(r["ref_norm"])) + return out + + +def ratio_vs(rows, key: str) -> dict[float, list[float]]: + """`{x: [monoprop/ppvm error ratio, ...]}` as `key` is swept. + + A ratio is only formed inside one fully-specified cell — same model, width, + depth, angle, seed and `atol` — so sweeping one axis never divides across + another. Remaining axes (seeds) pool into the list. + """ + cells: dict[tuple, dict[str, float]] = defaultdict(dict) + for r in rows: + if r["library"] not in PAIR: + continue + cell = (r["model"], r["qubits"], r["steps"], r["dt"], r["seed"], r["atol"]) + cells[cell][r["library"]] = float(r["l2_err"]) + cells[cell]["x"] = float(r[key]) + out: dict[float, list[float]] = defaultdict(list) + for libs in cells.values(): + if len(libs) < 3 or not libs.get(BASELINE): + continue + out[libs["x"]].append(libs["monoprop"] / libs[BASELINE]) + return dict(sorted(out.items())) + + +def band(ax, data, color, marker, label) -> None: + """Median line, with a min/max ribbon where an axis pooled several seeds.""" + xs = list(data) + ax.plot( + xs, + [statistics.median(data[x]) for x in xs], + color=color, + marker=marker, + markersize=5, + linewidth=1.6, + label=label, + ) + if any(len(v) > 1 for v in data.values()): + ax.fill_between( + xs, + [min(data[x]) for x in xs], + [max(data[x]) for x in xs], + color=color, + alpha=0.16, + linewidth=0, + ) + + +def panel_atol(ax, rows) -> None: + """Absolute error against atol, every engine, both Trotter models.""" + for model, dash in (("tfim", "-"), ("heisenberg", "--")): + for lib, (label, color, marker) in SERIES.items(): + pts = sorted( + (float(r["atol"]), float(r["l2_err"]) / float(r["ref_norm"])) + for r in rows + if r["model"] == model and r["library"] == lib + ) + if pts: + ax.plot( + *zip(*pts), + dash, + color=color, + marker=marker, + markersize=4, + linewidth=1.4, + alpha=0.9, + label=label if model == "tfim" else None, + ) + # Slope-1 guide, anchored on the baseline's own worst point so it sits beside + # the data rather than floating below it -- the eye is comparing gradients. + anchor = max( + (float(r["atol"]), float(r["l2_err"]) / float(r["ref_norm"])) + for r in rows + if r["model"] == "tfim" and r["library"] == BASELINE + ) + lo = min(float(r["atol"]) for r in rows) + ax.plot( + [lo, anchor[0]], + [anchor[1] * lo / anchor[0], anchor[1]], + ":", + color=TEXT_SECONDARY, + linewidth=1.3, + label="slope 1 (error ∝ atol)", + ) + ax.set_xscale("log") + ax.set_yscale("log") + ax.invert_xaxis() + ax.set_xlabel("truncation threshold atol") + ax.set_ylabel("relative L2 error ‖c − c*‖ / ‖c*‖") + ax.set_title( + "1 · Every rule converges at the same order\n" + "solid TFIM, dashed Heisenberg — n=8, 10 steps, dt=0.1", + fontsize=10, + color=TEXT_PRIMARY, + loc="left", + ) + ax.legend(fontsize=7, frameon=False, labelcolor=TEXT_SECONDARY, loc="lower left") + + +def panel_ratio(ax, rows, key, xlabel, title) -> None: + """monoprop / ppvm error ratio against one swept axis, one line per atol. + + The swept values are logarithmically spaced and uneven, so the axis is log + and ticks are pinned to the values actually run. + """ + ax.set_xscale("log") + ax.minorticks_off() + # Loosest threshold first: that is the series with the largest effect, and + # `sorted` on the raw strings would put 1e-5 ahead of 1e-3. + xs: list[float] = [] + for i, atol in enumerate( + sorted({r["atol"] for r in rows}, key=float, reverse=True) + ): + data = ratio_vs([r for r in rows if r["atol"] == atol], key) + if data: + band( + ax, + data, + RATIO_INK[i % len(RATIO_INK)], + RATIO_MARKS[i % len(RATIO_MARKS)], + f"atol = {atol}", + ) + xs = list(data) + # Set after the scale, which would otherwise install its own locator. + ax.set_xticks(xs) + ax.set_xticklabels([f"{x:g}" for x in xs], fontsize=8) + ax.axhline(1.0, color=SERIES["ppvm"][1], linewidth=1.2, linestyle=":") + ax.annotate( + "equal accuracy", + (0.03, 1.0), + xycoords=("axes fraction", "data"), + fontsize=7, + color=SERIES["ppvm"][1], + va="bottom", + ) + ax.set_xlabel(xlabel) + ax.set_ylabel("monoprop error / ppvm error") + ax.set_title(title, fontsize=10, color=TEXT_PRIMARY, loc="left") + ax.legend(fontsize=7, frameon=False, labelcolor=TEXT_SECONDARY) + + +def panel_scramble(ax, rows) -> None: + """The scrambling instance, small angle against the theta = pi/4 control.""" + for dt, dash, tag in ( + ("0.05", "-", "dt=0.05, 3200 gates"), + ("0.7853981633974483", "--", "dt=π/4, 80 gates"), + ): + sub = [r for r in rows if r["dt"] == dt] + for lib in PAIR: + data = relative(sub, lib) + if data: + label, color, marker = SERIES[lib] + ax.plot( + sorted(data), + [statistics.median(data[x]) for x in sorted(data)], + dash, + color=color, + marker=marker, + markersize=5, + linewidth=1.6, + label=f"{label} — {tag}", + ) + ax.set_xscale("log") + ax.set_yscale("log") + ax.invert_xaxis() + ax.set_xlabel("truncation threshold atol") + ax.set_ylabel("relative L2 error") + ax.set_title( + "4 · Scrambling instance: where they separate\n" + "n=8, all 65535 Pauli words populated, median of 5 seeds", + fontsize=10, + color=TEXT_PRIMARY, + loc="left", + ) + ax.legend(fontsize=7, frameon=False, labelcolor=TEXT_SECONDARY, loc="lower left") + + +CAPTION = ( + "Error is the L2 distance of the whole propagated coefficient vector from a " + "converged reference (ppvm at atol=1e-16), relative to ‖c*‖. ppvm, " + "PauliPropagation.jl and PauliStrings.jl accumulate every contribution to a term " + "and then threshold the merged coefficient; pauli-prop and monoprop threshold each " + "branch as it is emitted. A Pauli rotation sends Q to exactly two words, so a " + "merged value can exceed the larger contribution by at most 2x — which bounds the " + "divergence and ties both errors to the same atol. Panels 2-4: monoprop needs small " + "sin θ (a wide rejection band), depth (to populate and compound it) and a scrambled " + "operator (so the loss does not average out) simultaneously before it separates; " + "drop any one and the gap collapses or reverses. Ribbons are min/max over seeds." +) + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("--csv", type=Path, required=True, help="the 5-engine atol sweep") + ap.add_argument( + "--divergence", + type=Path, + required=True, + help="the depth / angle / scramble sweeps", + ) + ap.add_argument("--out", type=Path, required=True) + args = ap.parse_args() + + div = load(args.divergence) + # The two Heisenberg sweeps are told apart by width, which is how they were + # run: depth at n=6 (cheap enough for 320 steps), angle at n=8. + depth = [r for r in div if r["model"] == "heisenberg" and r["qubits"] == "6"] + angle = [r for r in div if r["model"] == "heisenberg" and r["qubits"] == "8"] + scramble = [r for r in div if r["model"] == "scramble"] + for name, rows in (("depth", depth), ("angle", angle), ("scramble", scramble)): + if not rows: + raise SystemExit(f"no {name} rows in {args.divergence}") + + fig, axes = plt.subplots(2, 2, figsize=(13.0, 9.6), facecolor=SURFACE) + for ax in axes.flat: + style(ax) + + panel_atol(axes[0][0], load(args.csv)) + panel_ratio( + axes[0][1], + depth, + "steps", + "Trotter steps", + "2 · The gap compounds with depth\nHeisenberg n=6, dt=0.05 (sin θ ≈ 0.1)", + ) + panel_ratio( + axes[1][0], + angle, + "dt", + "dt (θ = 2·dt, so dt = π/8 gives sin θ = 0.707)", + "3 · …and reverses at large angles\nHeisenberg n=8, 10 steps", + ) + panel_scramble(axes[1][1], scramble) + + fig.text( + 0.008, + 0.012, + textwrap.fill(CAPTION, 168), + fontsize=7.5, + color=TEXT_SECONDARY, + va="bottom", + ) + fig.tight_layout(rect=(0, 0.072, 1, 1)) + args.out.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(args.out, dpi=200, facecolor=SURFACE) + print(f"wrote {args.out}") + + +if __name__ == "__main__": + main() diff --git a/benchmarks/cross-library/plot_xbench.py b/benchmarks/cross-library/plot_xbench.py new file mode 100644 index 000000000..cd5f28dd5 --- /dev/null +++ b/benchmarks/cross-library/plot_xbench.py @@ -0,0 +1,244 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 +"""Render the cross-library benchmark from `run_xbench.py`'s `results.csv`. + +One row per model, three panels: runtime vs qubit count (log-y), runtime +relative to `ppvm`, and the size of the propagated operator — which is the +workload, and the thing that makes the runtimes comparable in the first place. + + uv run --no-project --with matplotlib python3 plot_xbench.py \ + --csv target/xbench/results.csv --out target/xbench/xbench.png +""" + +from __future__ import annotations + +import argparse +import csv +import textwrap +from collections import defaultdict +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +# Slots from the validated default palette, in this fixed order — it clears +# every hard gate on the adjacent pairlist that line charts use (worst CVD +# ΔE 9.2, worst normal-vision ΔE 27.5). Colour follows the entity, so a run with +# a subset of `--libs` keeps every survivor's hue. Magenta and aqua sit below +# 3:1 on the light surface, which obligates relief — `results.csv`, the driver's +# summary table, and the per-series markers are it. +BASELINE = "ppvm" +SERIES = { + "ppvm": ("ppvm (this repo)", "#2a78d6", "o"), + "pauli-propagation-jl": ("PauliPropagation.jl", "#eb6834", "s"), + "pauli-strings-jl": ("PauliStrings.jl", "#1baf7a", "^"), + "pauli-prop": ("pauli-prop (Qiskit)", "#4a3aa7", "D"), + "monoprop": ("monoprop", "#e87ba4", "v"), +} +TEXT_PRIMARY = "#0b0b0b" +TEXT_SECONDARY = "#52514e" +GRID = "#d8d7d3" +SURFACE = "#fcfcfb" +MODEL_TITLES = { + "tfim": "TFIM Trotter — ⟨0|Σ Z_i(t)|0⟩", + "heisenberg": "Heisenberg correlations — tr[Z₀ Z₀(t)]/2ⁿ", +} + + +def load(path: Path): + """`{model: {library: {n: row}}}`.""" + out: dict[str, dict[str, dict[int, dict[str, str]]]] = defaultdict( + lambda: defaultdict(dict) + ) + with path.open() as fh: + for row in csv.DictReader(fh): + out[row["model"]][row["library"]][int(row["qubits"])] = row + return out + + +def style(ax) -> None: + """Recessive grid and axes; the data carries the ink.""" + ax.set_facecolor(SURFACE) + ax.grid(True, which="major", color=GRID, linewidth=0.8, alpha=0.9) + ax.grid(True, which="minor", color=GRID, linewidth=0.5, alpha=0.5) + ax.set_axisbelow(True) + for side in ("top", "right"): + ax.spines[side].set_visible(False) + for side in ("left", "bottom"): + ax.spines[side].set_color(GRID) + ax.tick_params(colors=TEXT_SECONDARY, labelsize=9) + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("--csv", type=Path, required=True) + ap.add_argument("--out", type=Path, required=True) + ap.add_argument("--title", default=None) + args = ap.parse_args() + + data = load(args.csv) + models = [m for m in ("tfim", "heisenberg") if m in data] + fig, axes = plt.subplots( + len(models), + 3, + figsize=(15.5, 4.4 * len(models)), + squeeze=False, + facecolor=SURFACE, + ) + + for r, model in enumerate(models): + per_lib = data[model] + ax_t, ax_s, ax_n = axes[r] + widths = sorted({n for lib in per_lib.values() for n in lib}) + term_series: dict[str, tuple[int, ...]] = {} + for key, (label, color, marker) in SERIES.items(): + if key not in per_lib: + continue + ns = sorted(per_lib[key]) + times = [float(per_lib[key][n]["time_s"]) for n in ns] + ax_t.plot( + ns, + times, + color=color, + marker=marker, + markersize=6, + linewidth=2, + label=label, + ) + terms = [int(per_lib[key][n]["terms"]) for n in ns] + ax_n.plot( + ns, + terms, + color=color, + marker=marker, + markersize=6, + linewidth=2, + label=label, + ) + term_series[key] = tuple(terms) + if key != BASELINE and BASELINE in per_lib: + shared = [n for n in ns if n in per_lib[BASELINE]] + ratios = [ + float(per_lib[key][n]["time_s"]) + / float(per_lib[BASELINE][n]["time_s"]) + for n in shared + ] + ax_s.plot( + shared, + ratios, + color=color, + marker=marker, + markersize=6, + linewidth=2, + label=label, + ) + + # Engines whose support is identical draw the same line, so the ones + # underneath are invisible. Say so rather than let the panel imply that + # only the top series was measured. + if BASELINE in term_series: + same = [k for k, v in term_series.items() if v == term_series[BASELINE]] + if len(same) > 1: + ax_n.annotate( + "identical (exactly):\n" + + "\n".join(SERIES[k][0].split(" (")[0] for k in same), + # The support rises left-to-right, so the low-right corner + # is the one that stays clear of the marks. + xy=(0.97, 0.04), + xycoords="axes fraction", + fontsize=8, + color=TEXT_SECONDARY, + ha="right", + va="bottom", + ) + + ax_s.axhline(1.0, color=SERIES[BASELINE][1], linewidth=2, linestyle=(0, (4, 3))) + ax_s.annotate( + f"{SERIES[BASELINE][0]} = 1", + xy=(0.02, 1.0), + xycoords=("axes fraction", "data"), + va="bottom", + fontsize=8, + color=TEXT_SECONDARY, + ) + + for ax, ylabel, title in ( + (ax_t, "runtime (s, min of repeats)", "runtime"), + (ax_s, f"× {BASELINE}", "relative runtime"), + (ax_n, "terms in the propagated operator", "workload size"), + ): + style(ax) + ax.set_yscale("log") + # Qubit counts are integers; let matplotlib interpolate ticks and it + # invents 6.25 qubits. + ax.set_xticks(widths) + ax.set_xticklabels([str(n) for n in widths]) + ax.set_xlabel("qubits", color=TEXT_SECONDARY, fontsize=9) + ax.set_ylabel(ylabel, color=TEXT_SECONDARY, fontsize=9) + ax.set_title(title, color=TEXT_PRIMARY, fontsize=10, loc="left") + + ax_t.set_ylabel("runtime (s, min of repeats)", color=TEXT_SECONDARY, fontsize=9) + ax_t.annotate( + MODEL_TITLES.get(model, model), + xy=(0, 1.16), + xycoords="axes fraction", + fontsize=12, + fontweight="bold", + color=TEXT_PRIMARY, + ) + + # One figure-level legend above the panels: repeating it per row wastes + # space, and inside the axes it lands on the fastest series. + handles, labels = axes[0][0].get_legend_handles_labels() + fig.legend( + handles, + labels, + frameon=False, + fontsize=9.5, + labelcolor=TEXT_SECONDARY, + ncol=len(labels), + loc="upper left", + bbox_to_anchor=(0.006, 0.995), + ) + + sample = next(iter(next(iter(data.values())).values())) + row = next(iter(sample.values())) + n_engines = len({lib for per_lib in data.values() for lib in per_lib}) + caption = ( + f"first-order Trotter, steps={row['steps']}, dt={row['dt']}, " + f"truncation |c| < {row['atol']}; single-threaded; min of repeats. " + f"All {n_engines} engines are validated to propagate the identical operator " + "term-for-term before timing, but pauli-prop and monoprop prune at " + "branch-creation time rather than after accumulation, so at this " + "truncation they carry a different support — compare the workload panel " + "before reading their ratios. monoprop is capped to one thread " + "(monoprop_NUM_THREADS=1); uncapped it takes one partition per core. " + "Per-point numbers in results.csv." + ) + if args.title: + fig.suptitle( + args.title, fontsize=14, color=TEXT_PRIMARY, x=0.006, y=0.995, ha="left" + ) + # Legend sits just under the title when there is one. + fig.legends[0].set_bbox_to_anchor((0.006, 0.962)) + wrapped = textwrap.fill(caption, width=178) + fig.text( + 0.008, + 0.004, + wrapped, + fontsize=8.5, + color=TEXT_SECONDARY, + ha="left", + va="bottom", + ) + top = 0.90 if args.title else 0.945 + fig.tight_layout(rect=(0, 0.018 * (wrapped.count("\n") + 1) + 0.012, 1, top)) + args.out.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(args.out, dpi=160, facecolor=SURFACE) + print(f"wrote {args.out}") + + +if __name__ == "__main__": + main() diff --git a/benchmarks/cross-library/run_xbench.py b/benchmarks/cross-library/run_xbench.py new file mode 100644 index 000000000..cf1b9c5b0 --- /dev/null +++ b/benchmarks/cross-library/run_xbench.py @@ -0,0 +1,366 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 +"""Drive the cross-library Pauli-propagation benchmark and merge the results. + +Runs the same two circuits — TFIM Trotter and Heisenberg-model correlations — +through `ppvm`, PauliPropagation.jl, PauliStrings.jl, Qiskit's `pauli-prop` and +monoprop, all from one parameter contract, and writes one tidy CSV. + +Before it will report a single timing it **validates**: every engine dumps its +propagated support at a small width and the driver diffs them term-for-term +against `ppvm`. A cross-library benchmark whose engines are quietly computing +different things is worse than no benchmark, and every bug this check caught +while the harness was being written produced plausible-looking numbers that +agreed on the observable to nine digits — a reversed circuit in the `pauli-prop` +runner, duplicate Paulis collapsing in its readout, and monoprop's `exp(+iθH)` +sign convention. + + uv run --no-project python3 run_xbench.py --help +""" + +from __future__ import annotations + +import argparse +import csv +import os +import shutil +import subprocess +import sys +from dataclasses import dataclass, field +from pathlib import Path + +REPO = Path(__file__).resolve().parents[2] +JULIA_PROJECT = REPO / "julia-benchmarks" +BASELINE = "ppvm" +CSV_COLUMNS = [ + "model", + "library", + "qubits", + "steps", + "dt", + "atol", + "time_s", + "terms", + "observable", +] + + +@dataclass(frozen=True) +class Runner: + """One engine: how to invoke it, and whether it is available here.""" + + name: str + argv: list[str] + needs: str # the executable that must exist on PATH + env: dict[str, str] = field(default_factory=dict) + + def available(self) -> bool: + return shutil.which(self.needs) is not None + + +def _here(name: str) -> str: + return str(Path(__file__).with_name(name)) + + +RUNNERS = { + "ppvm": Runner( + "ppvm", + [ + "cargo", + "run", + "--release", + "-q", + "-p", + "ppvm-pauli-sum", + "--example", + "xbench", + ], + "cargo", + ), + "pauli-propagation-jl": Runner( + "pauli-propagation-jl", + [ + "julia", + f"--project={JULIA_PROJECT}", + "-t1", + str(JULIA_PROJECT / "benches" / "xbench_pp.jl"), + ], + "julia", + ), + "pauli-strings-jl": Runner( + "pauli-strings-jl", + [ + "julia", + f"--project={JULIA_PROJECT}", + "-t1", + str(JULIA_PROJECT / "benches" / "xbench_ps.jl"), + ], + "julia", + ), + "pauli-prop": Runner( + "pauli-prop", + [ + "uv", + "run", + "--no-project", + "--with", + "pauli-prop", + "python3", + _here("xbench_qiskit.py"), + ], + "uv", + ), + # monoprop takes one serial partition per physical core when left alone, so + # the cap is part of the invocation, not an optional extra. The runner + # re-asserts it from its own CPU/wall ratio and fails if it did not land. + "monoprop": Runner( + "monoprop", + [ + "uv", + "run", + "--no-project", + "--with", + "monoprop", + "python3", + _here("xbench_monoprop.py"), + ], + "uv", + {"monoprop_NUM_THREADS": "1", "monoprop_PARTITIONS": "off"}, + ), +} + + +def invoke(runner: Runner, env_extra: dict[str, str], quiet: bool) -> str: + """Run one engine and return its stdout, streaming its stderr as progress.""" + env = os.environ.copy() + env.update(env_extra) + env.update(runner.env) + proc = subprocess.run( + runner.argv, + cwd=REPO, + env=env, + capture_output=True, + text=True, + check=False, + ) + if proc.returncode != 0: + sys.stderr.write(proc.stderr) + raise SystemExit(f"{runner.name} failed with exit code {proc.returncode}") + if not quiet and proc.stderr.strip(): + sys.stderr.write(proc.stderr) + return proc.stdout + + +def parse_dump(text: str) -> dict[str, float]: + """Parse a `DUMP=1` support listing into `{word: coefficient}`.""" + terms: dict[str, float] = {} + for line in text.splitlines(): + if not line or line.startswith("#"): + continue + word, coeff = line.split() + terms[word] = float(coeff) + return terms + + +def validate( + libs: list[str], params: dict[str, str], atol_cmp: float, quiet: bool +) -> None: + """Assert every engine propagates the *same* operator, term for term.""" + print("validating: all engines must agree term-for-term", file=sys.stderr) + for model in ("tfim", "heisenberg"): + env = dict(params) + env.update( + {"MODEL": model, "QUBITS": "4", "STEPS": "3", "ATOL": "1e-14", "DUMP": "1"} + ) + reference: dict[str, float] | None = None + ref_lib = "" + for lib in libs: + terms = parse_dump(invoke(RUNNERS[lib], env, quiet=True)) + if reference is None: + reference, ref_lib = terms, lib + print( + f" {model}: {lib} — {len(terms)} terms (reference)", + file=sys.stderr, + ) + continue + missing = set(reference) - set(terms) + extra = set(terms) - set(reference) + worst = max( + (abs(reference[w] - terms[w]) for w in set(reference) & set(terms)), + default=0.0, + ) + if missing or extra or worst > atol_cmp: + raise SystemExit( + f"{model}: {lib} disagrees with {ref_lib} — " + f"{len(terms)} vs {len(reference)} terms, " + f"{len(missing)} missing, {len(extra)} extra, max|Δc|={worst:.3e} " + f"(tolerance {atol_cmp:.0e}). Refusing to report timings." + ) + print( + f" {model}: {lib} — {len(terms)} terms, max|Δc|={worst:.1e} OK", + file=sys.stderr, + ) + print("validation passed\n", file=sys.stderr) + + +def main() -> None: + ap = argparse.ArgumentParser( + description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter + ) + ap.add_argument("--model", default="both", choices=["tfim", "heisenberg", "both"]) + ap.add_argument( + "--qubits", default="8,12,16,20,24,28,32", help="applies to every model" + ) + ap.add_argument( + "--qubits-tfim", + default=None, + help="override --qubits for the TFIM sweep (its support grows linearly in n, " + "so it reaches far wider systems than the Heisenberg one at equal cost)", + ) + ap.add_argument( + "--qubits-heisenberg", default=None, help="override --qubits for Heisenberg" + ) + ap.add_argument("--steps", type=int, default=10) + ap.add_argument("--dt", type=float, default=0.1) + ap.add_argument("--j", type=float, default=1.0, help="bond coupling J") + ap.add_argument("--h", type=float, default=1.0, help="field strength h") + ap.add_argument( + "--atol", type=float, default=1e-6, help="coefficient truncation threshold" + ) + ap.add_argument( + "--iters", type=int, default=3, help="timed repeats; the minimum is reported" + ) + ap.add_argument("--libs", default=",".join(RUNNERS)) + ap.add_argument("--out", default="target/xbench", type=Path) + ap.add_argument("--skip-validate", action="store_true") + ap.add_argument( + "--validate-tol", + type=float, + default=1e-10, + help="absolute coefficient tolerance for the cross-engine agreement check", + ) + ap.add_argument( + "--reuse", + type=Path, + default=None, + help="re-print the summary from an existing results.csv", + ) + ap.add_argument("-q", "--quiet", action="store_true") + args = ap.parse_args() + + if args.reuse is not None: + with args.reuse.open() as fh: + summarize(list(csv.DictReader(fh))) + return + + libs = [lib.strip() for lib in args.libs.split(",") if lib.strip()] + for lib in libs: + if lib not in RUNNERS: + raise SystemExit(f"unknown library {lib!r}; known: {', '.join(RUNNERS)}") + missing = [lib for lib in libs if not RUNNERS[lib].available()] + if missing: + print( + f"skipping (missing {', '.join(RUNNERS[m].needs for m in missing)}): " + f"{', '.join(missing)}", + file=sys.stderr, + ) + libs = [lib for lib in libs if lib not in missing] + if not libs: + raise SystemExit("no runnable engines") + + params = { + "DT": repr(args.dt), + "JCOUP": repr(args.j), + "HFIELD": repr(args.h), + "ITERS": str(args.iters), + } + if not args.skip_validate: + validate(libs, params, args.validate_tol, args.quiet) + + models = ["tfim", "heisenberg"] if args.model == "both" else [args.model] + per_model_qubits = { + "tfim": args.qubits_tfim or args.qubits, + "heisenberg": args.qubits_heisenberg or args.qubits, + } + rows: list[dict[str, str]] = [] + for model in models: + for lib in libs: + env = dict(params) + env.update( + { + "MODEL": model, + "QUBITS": per_model_qubits[model], + "STEPS": str(args.steps), + "ATOL": repr(args.atol), + } + ) + print(f"running {model} / {lib}", file=sys.stderr) + out = invoke(RUNNERS[lib], env, args.quiet) + reader = csv.DictReader(line for line in out.splitlines() if line.strip()) + rows.extend(dict(row) for row in reader) + + out_dir = REPO / args.out if not args.out.is_absolute() else args.out + out_dir.mkdir(parents=True, exist_ok=True) + csv_path = out_dir / "results.csv" + with csv_path.open("w", newline="") as fh: + writer = csv.DictWriter(fh, fieldnames=CSV_COLUMNS) + writer.writeheader() + writer.writerows(rows) + + summarize(rows) + print(f"\nwrote {csv_path}", file=sys.stderr) + + +def summarize(rows: list[dict[str, str]]) -> None: + """Print a time table plus the speedup of every engine relative to `ppvm`.""" + by_model: dict[str, dict[int, dict[str, dict[str, str]]]] = {} + for row in rows: + by_model.setdefault(row["model"], {}).setdefault(int(row["qubits"]), {})[ + row["library"] + ] = row + + for model, per_n in by_model.items(): + libs = sorted({lib for cells in per_n.values() for lib in cells}) + print(f"\n=== {model} ===") + header = f"{'n':>4} {'terms':>9} " + " ".join(f"{lib:>22}" for lib in libs) + print(header) + print("-" * len(header)) + for n in sorted(per_n): + cells = per_n[n] + base = cells.get(BASELINE) + terms = base["terms"] if base else next(iter(cells.values()))["terms"] + line = f"{n:>4} {terms:>9} " + parts = [] + for lib in libs: + row = cells.get(lib) + if row is None: + parts.append(f"{'—':>22}") + continue + t = float(row["time_s"]) + if base is not None and lib != BASELINE: + ratio = t / float(base["time_s"]) + parts.append(f"{t:>11.4f}s ({ratio:5.2f}x)") + else: + parts.append(f"{t:>11.4f}s {'(1.00x)':>9}") + print(line + " ".join(parts)) + + # A runtime ratio only means something at equal work. Any engine whose + # truncation rule leaves it carrying a materially different support is + # called out here rather than left to the reader to spot in the plot. + drift: dict[str, list[str]] = {} + for n in sorted(per_n): + base = per_n[n].get(BASELINE) + if base is None: + continue + for lib, row in per_n[n].items(): + if lib == BASELINE: + continue + rel = int(row["terms"]) / int(base["terms"]) - 1.0 + if abs(rel) > 0.02: + drift.setdefault(lib, []).append(f"n={n}: {rel:+.0%}") + for lib, notes in drift.items(): + print(f" note: {lib} support vs {BASELINE} — {', '.join(notes)}") + + +if __name__ == "__main__": + main() diff --git a/benchmarks/cross-library/xbench_accuracy.py b/benchmarks/cross-library/xbench_accuracy.py new file mode 100644 index 000000000..4132936e2 --- /dev/null +++ b/benchmarks/cross-library/xbench_accuracy.py @@ -0,0 +1,197 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 +"""Measure what each engine's truncation rule costs in accuracy, not in time. + +`run_xbench.py` answers "how fast", holding `atol` fixed. It cannot answer "how +wrong", because the engines do not all mean the same thing by `atol`: `ppvm`, +PauliPropagation.jl and PauliStrings.jl accumulate every contribution to a term +and then drop the merged coefficient, while `pauli-prop` and monoprop test a +branch's prospective coefficient before they emit it. Same threshold, different +rule, so the fair question is how far each lands from the untruncated answer. + +This sweeps `atol` at a small width and diffs each engine's **whole coefficient +vector** against a converged reference. The vector norm is the point. Judging +this on the scalar `observable` is actively misleading, because its truncation +error changes sign — `ppvm` on TFIM walks 7.7e-4, 1.8e-3, 1.0e-4, 3.0e-7, 1.2e-6 +over `atol` 1e-3 … 1e-7, and that accidental near-zero at 1e-6 manufactures a 27x +gap against monoprop where the norm shows 6 %. Both are recorded so the +discrepancy stays visible rather than being something you have to rediscover. + +The reference is the baseline engine at `--ref-atol`. Widths whose reachable +sector saturates give an exact reference (Heisenberg `n=8` tops out at 16 384 +terms); elsewhere check it has converged before trusting a small ratio. + + uv run --no-project python3 xbench_accuracy.py --help +""" + +from __future__ import annotations + +import argparse +import csv +import itertools +import math +import sys +from pathlib import Path + +from run_xbench import BASELINE, RUNNERS, invoke, parse_dump + +CSV_COLUMNS = [ + "model", + "library", + "qubits", + "steps", + "dt", + "seed", + "atol", + "ref_atol", + "ref_terms", + "ref_norm", + "terms", + "l2_err", + "l1_err", + "max_coeff_err", + "lost_mass", + "kept_subthreshold", + "dropped_above_atol", + "observable", + "ref_observable", + "obs_abs_err", +] + + +def norm(terms: dict[str, float]) -> float: + """`‖c‖₂`. Unitary conjugation preserves it, so it is the scale the absolute + errors should be read against.""" + return math.sqrt(sum(v * v for v in terms.values())) + + +def measure(approx: dict[str, float], ref: dict[str, float], atol: float) -> dict: + """Compare one propagated support against the converged one.""" + words = set(approx) | set(ref) + diffs = [abs(approx.get(w, 0.0) - ref.get(w, 0.0)) for w in words] + return { + "terms": len(approx), + "l2_err": math.sqrt(sum(d * d for d in diffs)), + "l1_err": sum(diffs), + "max_coeff_err": max(diffs, default=0.0), + # Weight the engine discarded that the reference says was there. + "lost_mass": sum(abs(c) for w, c in ref.items() if w not in approx), + # Rows held whose converged value is under the threshold they were given: + # the signature of a branch that cleared `atol` when it was emitted and + # was never re-tested after later contributions cancelled it. + "kept_subthreshold": sum(1 for w in approx if abs(ref.get(w, 0.0)) < atol), + # The opposite failure: mass thrown away that belonged above `atol`. + "dropped_above_atol": sum( + 1 for w, c in ref.items() if w not in approx and abs(c) >= atol + ), + } + + +def observable(text: str) -> float: + """Pull the `observable` column out of a runner's one-row CSV.""" + return float( + next(csv.DictReader(l for l in text.splitlines() if l.strip()))["observable"] + ) + + +def run(lib: str, env: dict[str, str], quiet: bool) -> tuple[dict[str, float], float]: + """Get one engine's support and its readout at these parameters.""" + dump = parse_dump(invoke(RUNNERS[lib], {**env, "DUMP": "1"}, quiet)) + return dump, observable(invoke(RUNNERS[lib], env, quiet)) + + +def main() -> None: + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("--models", default="tfim,heisenberg") + ap.add_argument("--qubits", type=int, default=8) + ap.add_argument( + "--steps", + default="10", + help="comma-separated depths; the reference is recomputed per depth", + ) + ap.add_argument("--dt", default="0.1", help="comma-separated angles scales") + ap.add_argument("--atols", default="1e-3,1e-4,1e-5,1e-6,1e-7") + ap.add_argument( + "--seeds", + default="12345", + help="comma-separated circuit seeds; only `scramble` reads them, and it " + "needs several because a random instance varies", + ) + ap.add_argument( + "--ref-atol", + default="1e-16", + help="threshold for the reference run; must be converged at this width", + ) + ap.add_argument("--libs", default=",".join(RUNNERS)) + ap.add_argument("--out", type=Path, default=Path("target/xbench/accuracy.csv")) + ap.add_argument("--quiet", action="store_true") + args = ap.parse_args() + + libs = [lib for lib in args.libs.split(",") if lib.strip()] + if unknown := set(libs) - set(RUNNERS): + raise SystemExit(f"unknown libraries: {sorted(unknown)}") + if missing := [lib for lib in libs if not RUNNERS[lib].available()]: + print(f"skipping (toolchain absent): {', '.join(missing)}", file=sys.stderr) + libs = [lib for lib in libs if lib not in missing] + if BASELINE not in libs: + raise SystemExit(f"{BASELINE} is the reference and cannot be skipped") + + base = {"QUBITS": str(args.qubits), "ITERS": "1", "JCOUP": "1.0", "HFIELD": "1.0"} + grid = itertools.product( + args.models.split(","), + args.steps.split(","), + args.dt.split(","), + args.seeds.split(","), + ) + rows = [] + for model, steps, dt, seed in grid: + env = {**base, "MODEL": model, "STEPS": steps, "DT": dt, "SEED": seed} + # The reference depends on every axis, so it is recomputed per cell. + ref, ref_obs = run(BASELINE, {**env, "ATOL": args.ref_atol}, args.quiet) + print( + f"{model} n={args.qubits} steps={steps} dt={dt} seed={seed}: reference is" + f" {BASELINE} at atol={args.ref_atol} — {len(ref)} terms," + f" observable {ref_obs:.12g}", + file=sys.stderr, + ) + for atol in args.atols.split(","): + for lib in libs: + approx, obs = run(lib, {**env, "ATOL": atol}, args.quiet) + stats = measure(approx, ref, float(atol)) + rows.append( + { + "model": model, + "library": lib, + "qubits": args.qubits, + "steps": steps, + "dt": dt, + "seed": seed, + "atol": atol, + "ref_atol": args.ref_atol, + "ref_terms": len(ref), + "ref_norm": f"{norm(ref):.6e}", + "observable": f"{obs:.15g}", + "ref_observable": f"{ref_obs:.15g}", + "obs_abs_err": f"{abs(obs - ref_obs):.6e}", + **{ + k: f"{v:.6e}" if isinstance(v, float) else v + for k, v in stats.items() + }, + } + ) + print( + f" atol={atol} {lib:>20}: {stats['terms']:>7} terms," + f" L2 {stats['l2_err']:.3e}, obs err {rows[-1]['obs_abs_err']}", + file=sys.stderr, + ) + + args.out.parent.mkdir(parents=True, exist_ok=True) + with args.out.open("w", newline="") as fh: + writer = csv.DictWriter(fh, fieldnames=CSV_COLUMNS) + writer.writeheader() + writer.writerows(rows) + print(f"wrote {len(rows)} rows to {args.out}", file=sys.stderr) + + +if __name__ == "__main__": + main() diff --git a/benchmarks/cross-library/xbench_monoprop.py b/benchmarks/cross-library/xbench_monoprop.py new file mode 100644 index 000000000..2a2c2abfc --- /dev/null +++ b/benchmarks/cross-library/xbench_monoprop.py @@ -0,0 +1,301 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 +"""The monoprop side of the cross-library Pauli-propagation benchmark. + +See `README.md` for the shared circuit definitions, the parameter contract, and +the CSV schema — every runner reads the same environment variables and prints the +same columns. + +Three things about monoprop need saying, because each of them silently changes +the number this file prints: + +**It is parallel by default.** With one MPI rank the engine takes one serial +partition per *physical core* (`resolve_partition_count_`), so an uncapped run on +this 14-core machine is a ~9.5x-CPU, ~3x-wall-faster run than the serial one, and +is not comparable with the other single-threaded engines. `monoprop_NUM_THREADS=1` +and `monoprop_PARTITIONS=off` pin it to one partition, and the runner **measures +its own CPU/wall ratio and fails** if the cap did not take. An environment +variable that is read once into a cached C++ static is exactly the kind of setting +that goes stale without anyone noticing. + +**Its angles are not the spec's angles.** `ExpGate` applies `exp(+iθH)`, so the +spec's `exp(-iθ/2·G)` needs `θ_monoprop = -θ_spec/2`. All four other engines take +`θ_spec` directly. + +**It reports more rows than it has terms.** The engine retains monomials whose +coefficient has cancelled to exactly zero, and `size()` counts those. At +`n=12` Heisenberg that is 3.06M tracked rows against 1.32M terms above the +truncation threshold. The `terms` column is the above-threshold support, which is +what the other four engines mean by it; the tracked size goes to stderr. + + MODEL=tfim QUBITS=8,16,24 STEPS=10 DT=0.1 JCOUP=1.0 HFIELD=1.0 ATOL=1e-6 \ + uv run --no-project --with monoprop python3 xbench_monoprop.py +""" + +from __future__ import annotations + +import os + +# Must precede the import: the C++ side reads these once into a cached static +# (`monoprop::config::get()`), so a later assignment would be ignored. The driver +# sets them in the subprocess environment too; this makes a direct invocation of +# this file single-threaded as well, rather than quietly parallel. +os.environ.setdefault("monoprop_NUM_THREADS", "1") +os.environ.setdefault("monoprop_PARTITIONS", "off") + +import resource +import sys +import time + +from monoprop import ( + Circuit, + ExpGate, + Pauli, + PauliOperator, + PauliPropagator, +) + +MODEL = os.environ.get("MODEL", "tfim") +STEPS = int(os.environ.get("STEPS", "10")) +DT = float(os.environ.get("DT", "0.1")) +JCOUP = float(os.environ.get("JCOUP", "1.0")) +HFIELD = float(os.environ.get("HFIELD", "1.0")) +ATOL = float(os.environ.get("ATOL", "1e-6")) +ITERS = int(os.environ.get("ITERS", "3")) +SEED = int(os.environ.get("SEED", "12345")) +# A serial run has CPU/wall ~1. Anything materially above it means the thread cap +# did not take and the timing is not comparable with the other engines. +CPU_WALL_MAX = float(os.environ.get("CPU_WALL_MAX", "1.5")) + +THETA_BOND = 2 * JCOUP * DT +THETA_SITE = 2 * HFIELD * DT + + +def cpu_seconds() -> float: + """Total CPU time (user+sys) charged to this process and its children.""" + me = resource.getrusage(resource.RUSAGE_SELF) + kids = resource.getrusage(resource.RUSAGE_CHILDREN) + return me.ru_utime + me.ru_stime + kids.ru_utime + kids.ru_stime + + +def seed_operator(n: int) -> PauliOperator: + """`Σ_i Z_i` for TFIM, `Z_0` for Heisenberg and scramble.""" + if MODEL == "tfim": + terms: dict[Pauli | str, float] = {Pauli("Z", (i,)): 1.0 for i in range(n)} + else: + terms = {Pauli("Z", (0,)): 1.0} + return PauliOperator(terms, n) + + +class SplitMix64: + """splitmix64, matching `SplitMix64` in `examples/xbench.rs` bit for bit. + + Both runners have to emit the *same* random circuit, and neither language's + stdlib RNG is specified tightly enough to rely on. This is short enough to + state in both and is checked by the term-for-term dump diff. + """ + + MASK = (1 << 64) - 1 + + def __init__(self, seed: int) -> None: + self.state = seed & self.MASK + + def next(self) -> int: + self.state = (self.state + 0x9E3779B97F4A7C15) & self.MASK + z = self.state + z = ((z ^ (z >> 30)) * 0xBF58476D1CE4E5B9) & self.MASK + z = ((z ^ (z >> 27)) * 0x94D049BB133111EB) & self.MASK + return z ^ (z >> 31) + + def unit(self) -> float: + """Uniform in `[0, 1)`, from the top 53 bits.""" + return (self.next() >> 11) / float(1 << 53) + + +AXES = "XYZ" + + +def scramble_gates(n: int) -> list[tuple[str, tuple[int, ...], float]]: + """`STEPS · n` random all-to-all two-qubit Pauli rotations. + + Draw order — pair, offset, axis, axis, angle — must match `scramble_gates` + in `examples/xbench.rs` exactly, including the offset trick that avoids a + rejection loop for `b != a` (a loop would consume a variable number of + draws and desynchronise the two streams). + """ + if n < 2: + raise SystemExit("scramble needs at least 2 qubits") + rng = SplitMix64(SEED) + theta_max = 2 * JCOUP * DT + gates = [] + for _ in range(STEPS * n): + a = rng.next() % n + b = (a + 1 + rng.next() % (n - 1)) % n + pa = AXES[rng.next() % 3] + pb = AXES[rng.next() % 3] + theta = theta_max * rng.unit() + # P_a ⊗ P_b is symmetric under swapping both, so normalise to ascending + # qubit order -- monoprop's `Pauli` expects sorted supports. + if a > b: + a, b, pa, pb = b, a, pb, pa + gates.append((pa + pb, (a, b), theta)) + return gates + + +def gate_list(n: int) -> list[tuple[str, tuple[int, ...], float]]: + """The gate sequence in *application* order, per the spec.""" + if MODEL == "scramble": + return scramble_gates(n) + gates: list[tuple[str, tuple[int, ...], float]] = [] + for _ in range(STEPS): + if MODEL == "tfim": + gates += [("X", (i,), THETA_SITE) for i in range(n)] + gates += [("ZZ", (i, i + 1), THETA_BOND) for i in range(n - 1)] + else: + for i in range(n - 1): + gates += [ + ("XX", (i, i + 1), THETA_BOND), + ("YY", (i, i + 1), THETA_BOND), + ("ZZ", (i, i + 1), THETA_BOND), + ] + gates += [("Z", (i,), THETA_SITE) for i in range(n)] + return gates + + +def build_circuit(n: int) -> Circuit: + """The spec's gate sequence as a monoprop `Circuit`. + + Two conversions happen here, and the term-for-term validation in + `run_xbench.py` is what pinned both of them down: + + * **Angle.** `ExpGate` applies `exp(+iθH)` (its own docstring flags the + positive sign as the difference from Qiskit's `r

`), while the shared + spec's gate is `exp(-iθ_spec/2 · G)`. Hence `θ = -θ_spec/2`. + * **Order.** Like Qiskit's `pauli-prop`, the propagator conjugates from the + end of the gate list backwards, so the spec's application order is + obtained by appending in reverse. + + Neither is cosmetic. All eight sign/order combinations were checked against + the reference dump: forward order loses terms outright (109 against TFIM's + 124), and every wrong angle keeps the right support while moving + coefficients by up to 1.3. Only this one lands within 5e-13. + """ + gates = list(reversed(gate_list(n))) + exp_gates = [ + ExpGate(PauliOperator({Pauli(string, qubits): 1.0}, n), index) + for index, (string, qubits, _theta) in enumerate(gates) + ] + parameters = [-0.5 * theta for (_s, _q, theta) in gates] + return Circuit(gates=exp_gates, parameters=parameters) + + +def propagate(n: int, circuit: Circuit) -> PauliPropagator: + """Heisenberg-picture propagation with coefficient-only truncation. + + `cutoff` is monoprop's mandatory bound on retained Pauli *weight*, which the + other four engines have no analogue of; `cutoff=n` is the whole register, so + it never binds and `lower_atol` is left as the only truncation — the shared + `|c| < atol` rule. + """ + return PauliPropagator.from_circuit( + circuit, seed_operator(n), cutoff=n, lower_atol=ATOL + ) + + +def support(mp: PauliPropagator, n: int) -> dict[str, float]: + """`{word: coefficient}` with site 0 leftmost, exact-zero rows dropped. + + `atol=0.0` asks the engine for everything it holds, and the threshold is + applied here so the filter is the same one the other engines apply during + propagation. + """ + operator = mp.evolved_operator(atol=0.0) + terms: dict[str, float] = {} + for pauli, coeff in operator.terms.items(): + value = float(coeff.real if isinstance(coeff, complex) else coeff) + if abs(value) < ATOL: + continue + chars = ["I"] * n + for qubit, letter in zip(pauli.qubits, pauli.string): + chars[qubit] = letter + terms["".join(chars)] = value + return terms + + +def readout(terms: dict[str, float], n: int) -> float: + """`⟨0…0|O|0…0⟩` for TFIM; the `Z_0` autocorrelator for Heisenberg.""" + if MODEL == "tfim": + # ⟨0|Z|0⟩ = 1 and ⟨0|X|0⟩ = ⟨0|Y|0⟩ = 0, so only the X/Y-free terms survive. + return sum(c for w, c in terms.items() if "X" not in w and "Y" not in w) + return terms.get("Z" + "I" * (n - 1), 0.0) + + +def dump(n: int) -> None: + """Print the propagated support as `word coefficient`, largest first.""" + terms = support(propagate(n, build_circuit(n)), n) + print(f"# {len(terms)} terms") + for word in sorted(terms, key=lambda w: (-abs(terms[w]), w)): + print(f"{word} {terms[word]:+.12e}") + + +def main() -> None: + qubits = [ + int(t) for t in os.environ.get("QUBITS", "8,12,16,20,24,28,32").split(",") + ] + if os.environ.get("DUMP"): + dump(qubits[0]) + return + print( + f"monoprop {MODEL}: steps={STEPS} dt={DT} J={JCOUP} h={HFIELD} atol={ATOL} " + f"iters={ITERS} threads={os.environ['monoprop_NUM_THREADS']} " + f"partitions={os.environ['monoprop_PARTITIONS']} cores={os.cpu_count()}", + file=sys.stderr, + ) + print("model,library,qubits,steps,dt,atol,time_s,terms,observable") + for n in qubits: + # Circuit construction is not propagation; built once, outside the timed + # region, as in every other runner. monoprop's `propagate` does expand + # the gate list internally each call, but that is 0.0–2.2% of the total + # at these widths (measured), not enough to move a ratio. + circuit = build_circuit(n) + + best = float("inf") + cpu_at_best = 0.0 + propagator = None + for _ in range(ITERS): + cpu0, wall0 = cpu_seconds(), time.perf_counter() + propagator = propagate(n, circuit) + wall = time.perf_counter() - wall0 + if wall < best: + best, cpu_at_best = wall, cpu_seconds() - cpu0 + + # A parallel run would report a wall time the other engines never had the + # chance to compete with. Refuse to print it. + ratio = cpu_at_best / best if best > 0 else 1.0 + if ratio > CPU_WALL_MAX: + raise SystemExit( + f"n={n}: monoprop used {ratio:.2f}x CPU per wall second " + f"({cpu_at_best:.3f}s CPU in {best:.3f}s wall), so the thread cap " + f"did not take and this timing is not comparable with the other " + f"single-threaded engines. Check monoprop_NUM_THREADS / " + f"monoprop_PARTITIONS reach the subprocess." + ) + + assert propagator is not None + # Materializing the support is O(size) in Python and dwarfs the + # propagation at the wide end, so it happens once, after timing. + terms = support(propagator, n) + obs = readout(terms, n) + print( + f"{MODEL},monoprop,{n},{STEPS},{DT},{ATOL},{best:.7g},{len(terms)},{obs!r}" + ) + print( + f" n={n} {best:.4f}s {len(terms)} terms obs={obs} " + f"[cpu/wall={ratio:.2f} tracked={propagator.size()}]", + file=sys.stderr, + ) + sys.stdout.flush() + + +if __name__ == "__main__": + main() diff --git a/benchmarks/cross-library/xbench_qiskit.py b/benchmarks/cross-library/xbench_qiskit.py new file mode 100644 index 000000000..6122a60e3 --- /dev/null +++ b/benchmarks/cross-library/xbench_qiskit.py @@ -0,0 +1,190 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 +"""The Qiskit `pauli-prop` side of the cross-library Pauli-propagation benchmark. + +See `README.md` for the shared circuit definitions, the parameter contract, and +the CSV schema — every runner reads the same environment variables and prints +the same columns. + +`pauli-prop` truncates on `atol` *and* on a mandatory `max_terms` cap, which the +other three engines do not have. `MAX_TERMS` therefore defaults high enough to +be non-binding, and the runner asserts it never bound: if the propagated support +ever reaches the cap the row is not comparable and the run fails rather than +quietly reporting a differently-truncated number. + + MODEL=tfim QUBITS=8,16,24 STEPS=10 DT=0.1 JCOUP=1.0 HFIELD=1.0 ATOL=1e-6 \ + uv run --no-project --with pauli-prop python3 xbench_qiskit.py +""" + +from __future__ import annotations + +import os +import sys +import time + +import numpy as np +from pauli_prop import propagate_through_circuit +from qiskit import QuantumCircuit +from qiskit.quantum_info import SparsePauliOp + +MODEL = os.environ.get("MODEL", "tfim") +STEPS = int(os.environ.get("STEPS", "10")) +DT = float(os.environ.get("DT", "0.1")) +JCOUP = float(os.environ.get("JCOUP", "1.0")) +HFIELD = float(os.environ.get("HFIELD", "1.0")) +ATOL = float(os.environ.get("ATOL", "1e-6")) +ITERS = int(os.environ.get("ITERS", "3")) +MAX_TERMS = int(os.environ.get("MAX_TERMS", str(1 << 22))) + +THETA_BOND = 2 * JCOUP * DT +THETA_SITE = 2 * HFIELD * DT + + +def seed_operator(n: int) -> SparsePauliOp: + """`Σ_i Z_i` for TFIM, `Z_0` for Heisenberg. + + Qiskit Pauli labels are little-endian: the rightmost character is qubit 0. + """ + if MODEL == "tfim": + labels = ["I" * (n - 1 - i) + "Z" + "I" * i for i in range(n)] + else: + labels = ["I" * (n - 1) + "Z"] + return SparsePauliOp(labels, coeffs=np.ones(len(labels))) + + +def gate_list(n: int) -> list[tuple[str, int]]: + """`STEPS` first-order Trotter steps in *application* order, per the spec.""" + gates: list[tuple[str, int]] = [] + for _ in range(STEPS): + if MODEL == "tfim": + gates += [("rx", i) for i in range(n)] + gates += [("rzz", i) for i in range(n - 1)] + else: + for i in range(n - 1): + gates += [("rxx", i), ("ryy", i), ("rzz", i)] + gates += [("rz", i) for i in range(n)] + return gates + + +def build_circuit(n: int) -> QuantumCircuit: + """The spec's gate sequence, as a circuit `pauli-prop` will apply in order. + + A Qiskit circuit whose instructions are appended `g_1 … g_k` denotes the + unitary `U = g_k ⋯ g_1`, and in the Heisenberg frame `pauli-prop` computes + `U† O U` by conjugating **from the end of the instruction list backwards**. + So the gate `pauli-prop` applies first is the one appended last, and the + spec's application order is obtained by appending in reverse. + + This is not cosmetic. Appending in forward order propagates a genuinely + different operator: at `n=4, steps=3, atol=1e-14` it yields 108 terms + against TFIM's reference 124 and 61 against Heisenberg's 64, with + coefficients off by up to 0.1. Reversed, all four engines agree + term-for-term to the last bit — which is what `run_xbench.py --validate` + checks before it will report a timing. + """ + qc = QuantumCircuit(n) + for name, i in reversed(gate_list(n)): + if name == "rx": + qc.rx(THETA_SITE, i) + elif name == "rz": + qc.rz(THETA_SITE, i) + elif name == "rxx": + qc.rxx(THETA_BOND, i, i + 1) + elif name == "ryy": + qc.ryy(THETA_BOND, i, i + 1) + elif name == "rzz": + qc.rzz(THETA_BOND, i, i + 1) + else: # pragma: no cover - the list above is closed + raise ValueError(f"unknown gate {name}") + return qc + + +def support(op: SparsePauliOp) -> dict[str, float]: + """`{word: coefficient}` with site 0 leftmost, duplicate Paulis summed. + + `propagate_through_circuit` returns an operator that may list the same Pauli + more than once — `len(op)` counts rows, not distinct terms — so the support + size comparable with the other three engines is the size of this mapping. + Qiskit labels are little-endian, hence the reversal. + """ + terms: dict[str, float] = {} + for pauli, coeff in zip(op.paulis, np.asarray(op.coeffs)): + word = str(pauli)[::-1] + terms[word] = terms.get(word, 0.0) + float(np.real(coeff)) + return terms + + +def readout(op: SparsePauliOp, n: int) -> float: + """`⟨0…0|O|0…0⟩` for TFIM; the `Z_0` autocorrelator for Heisenberg.""" + paulis = op.paulis + coeffs = np.asarray(op.coeffs) + if MODEL == "tfim": + # ⟨0|Z|0⟩ = 1 and ⟨0|X|0⟩ = ⟨0|Y|0⟩ = 0, so only the X-free terms survive. + diagonal = ~paulis.x.any(axis=1) + return float(np.real(coeffs[diagonal].sum())) + z0_x = np.zeros(n, dtype=bool) + z0_z = np.zeros(n, dtype=bool) + z0_z[0] = True + hit = (paulis.x == z0_x).all(axis=1) & (paulis.z == z0_z).all(axis=1) + return float(np.real(coeffs[hit].sum())) + + +def dump(n: int) -> None: + """Print the propagated support as `word coefficient`, largest first. + + Site 0 leftmost, matching the Rust runner's `DUMP=1` output so the driver + can diff the two term-for-term. + """ + out, _bias = propagate_through_circuit( + seed_operator(n), build_circuit(n), max_terms=MAX_TERMS, atol=ATOL, frame="h" + ) + terms = support(out) + print(f"# {len(terms)} terms") + for word in sorted(terms, key=lambda w: (-abs(terms[w]), w)): + print(f"{word} {terms[word]:+.12e}") + + +def main() -> None: + qubits = [ + int(t) for t in os.environ.get("QUBITS", "8,12,16,20,24,28,32").split(",") + ] + if os.environ.get("DUMP"): + dump(qubits[0]) + return + print( + f"pauli-prop {MODEL}: steps={STEPS} dt={DT} J={JCOUP} h={HFIELD} " + f"atol={ATOL} iters={ITERS} max_terms={MAX_TERMS}", + file=sys.stderr, + ) + print("model,library,qubits,steps,dt,atol,time_s,terms,observable") + for n in qubits: + # Circuit construction is not propagation; built once, outside the + # timed region, as in every other runner. + circuit = build_circuit(n) + seed = seed_operator(n) + + best = float("inf") + terms = 0 + obs = float("nan") + for _ in range(ITERS): + t0 = time.perf_counter() + out, _bias = propagate_through_circuit( + seed, circuit, max_terms=MAX_TERMS, atol=ATOL, frame="h" + ) + best = min(best, time.perf_counter() - t0) + # Distinct Paulis, not the returned operator's row count. + terms = len(support(out)) + obs = readout(out, n) + if terms >= MAX_TERMS: + raise SystemExit( + f"n={n}: support hit the max_terms cap ({terms} >= {MAX_TERMS}); " + "raise MAX_TERMS — this row would be truncated differently from " + "the other engines and is not comparable" + ) + print(f"{MODEL},pauli-prop,{n},{STEPS},{DT},{ATOL},{best:.7g},{terms},{obs!r}") + print(f" n={n} {best:.4f}s {terms} terms obs={obs}", file=sys.stderr) + sys.stdout.flush() + + +if __name__ == "__main__": + main() diff --git a/crates/ppvm-pauli-sum/examples/xbench.rs b/crates/ppvm-pauli-sum/examples/xbench.rs new file mode 100644 index 000000000..3b048487b --- /dev/null +++ b/crates/ppvm-pauli-sum/examples/xbench.rs @@ -0,0 +1,354 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! The `ppvm` side of the cross-library Pauli-propagation benchmark. +//! +//! Two Trotter workloads — TFIM magnetization and a Heisenberg autocorrelator — +//! propagated in the Heisenberg picture with a coefficient-magnitude truncation +//! after every gate. See `benchmarks/cross-library/README.md` for the shared +//! spec: the gate order, the `θ = 2·c·dt` convention, the environment contract, +//! and the CSV schema that every runner in that harness prints. +//! +//! ```bash +//! MODEL=tfim QUBITS=8,16,24 STEPS=10 DT=0.1 JCOUP=1.0 HFIELD=1.0 ATOL=1e-6 \ +//! cargo run --release -p ppvm-pauli-sum --example xbench +//! ``` + +use std::time::Instant; + +use ppvm_pauli_sum::prelude::*; +use ppvm_pauli_sum::strategy::CoefficientThreshold; + +#[derive(Clone, Copy)] +struct Params { + model: Model, + steps: usize, + dt: f64, + j: f64, + h: f64, + atol: f64, + iters: usize, + seed: u64, +} + +#[derive(Clone, Copy, PartialEq, Eq)] +enum Model { + Tfim, + Heisenberg, + Scramble, +} + +impl Model { + fn parse(s: &str) -> Self { + match s { + "tfim" => Model::Tfim, + "heisenberg" => Model::Heisenberg, + "scramble" => Model::Scramble, + other => { + panic!("unknown MODEL {other:?} (expected `tfim`, `heisenberg` or `scramble`)") + } + } + } + + fn name(self) -> &'static str { + match self { + Model::Tfim => "tfim", + Model::Heisenberg => "heisenberg", + Model::Scramble => "scramble", + } + } +} + +/// splitmix64. Reimplemented rather than pulled in as a dependency because the +/// monoprop runner has to emit a bit-identical gate sequence from Python, and +/// this is short enough to state twice and check against a term-for-term diff. +struct SplitMix64(u64); + +impl SplitMix64 { + fn next(&mut self) -> u64 { + self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15); + let mut z = self.0; + z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9); + z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB); + z ^ (z >> 31) + } + + /// Uniform in `[0, 1)`, from the top 53 bits. + fn unit(&mut self) -> f64 { + (self.next() >> 11) as f64 / (1u64 << 53) as f64 + } +} + +/// One two-qubit Pauli rotation `exp(-i θ/2 · P_a ⊗ P_b)`. +#[derive(Clone, Copy)] +struct Gate { + axis_a: [u8; 2], + axis_b: [u8; 2], + a: usize, + b: usize, + theta: f64, +} + +/// `[x, z]` bits for the three non-identity Paulis, indexed `0..3`. +const AXES: [[u8; 2]; 3] = [[1, 0], [1, 1], [0, 1]]; + +/// The `scramble` workload: `steps · n` two-qubit Pauli rotations on uniformly +/// random *all-to-all* pairs, with random axes and random angles in `(0, 2·J·dt]`. +/// +/// Unlike the two Trotter models this has no lattice, no conserved quantity and +/// no uniform angle, so the propagated operator spreads over the whole `4^n` +/// space and the coefficient distribution is genuinely scrambled rather than +/// hierarchically ordered by Pauli weight. +fn scramble_gates(n: usize, p: Params) -> Vec { + assert!(n >= 2, "scramble needs at least 2 qubits"); + let mut rng = SplitMix64(p.seed); + let theta_max = 2.0 * p.j * p.dt; + (0..p.steps * n) + .map(|_| { + let a = (rng.next() % n as u64) as usize; + // Offset by 1..n-1 so `b != a` without a rejection loop, which would + // desynchronise the two implementations' draw counts. + let b = (a + 1 + (rng.next() % (n as u64 - 1)) as usize) % n; + let axis_a = AXES[(rng.next() % 3) as usize]; + let axis_b = AXES[(rng.next() % 3) as usize]; + Gate { + axis_a, + axis_b, + a, + b, + theta: theta_max * rng.unit(), + } + }) + .collect() +} + +fn env_f64(key: &str, default: f64) -> f64 { + std::env::var(key) + .ok() + .and_then(|s| s.parse().ok()) + .unwrap_or(default) +} + +fn env_usize(key: &str, default: usize) -> usize { + std::env::var(key) + .ok() + .and_then(|s| s.parse().ok()) + .unwrap_or(default) +} + +impl Params { + fn from_env() -> Self { + Params { + model: Model::parse(&std::env::var("MODEL").unwrap_or_else(|_| "tfim".to_string())), + steps: env_usize("STEPS", 10), + dt: env_f64("DT", 0.1), + j: env_f64("JCOUP", 1.0), + h: env_f64("HFIELD", 1.0), + atol: env_f64("ATOL", 1e-6), + iters: env_usize("ITERS", 3), + seed: env_usize("SEED", 12345) as u64, + } + } +} + +/// Smallest power-of-two byte width that holds `n` qubits with room to spare. +fn storage_bytes(n: usize) -> usize { + let need = n.div_ceil(8); + let mut k = 0; + while (1usize << k) <= need { + k += 1; + } + 1usize << k +} + +/// A single-site Pauli word `p` at `site`, identity elsewhere. +fn site_word(n: usize, site: usize, p: char) -> String { + (0..n).map(|j| if j == site { p } else { 'I' }).collect() +} + +type Cfg = config::fxhash::Byte>; + +/// Build the seed observable for `model` on `n` sites. +fn seed(n: usize, p: Params) -> PauliSum> { + let mut sum: PauliSum> = PauliSum::builder() + .n_qubits(n) + .strategy(CoefficientThreshold(p.atol)) + .capacity(1 << 12) + .build(); + match p.model { + Model::Tfim => { + for i in 0..n { + sum += (PauliWord::from(site_word(n, i, 'Z').as_str()), 1.0); + } + } + Model::Heisenberg | Model::Scramble => { + sum += (PauliWord::from(site_word(n, 0, 'Z').as_str()), 1.0); + } + } + sum +} + +/// Propagate the shared gate sequence through `state`. +fn propagate(state: &mut PauliSum>, n: usize, p: Params, gates: &[Gate]) { + let theta_bond = 2.0 * p.j * p.dt; + let theta_site = 2.0 * p.h * p.dt; + if p.model == Model::Scramble { + for g in gates { + state.rotate_2(g.axis_a, g.axis_b, g.a, g.b, g.theta); + state.truncate(); + } + return; + } + for _ in 0..p.steps { + match p.model { + Model::Tfim => { + for i in 0..n { + state.rx(i, theta_site); + state.truncate(); + } + for i in 0..n.saturating_sub(1) { + state.rzz(i, i + 1, theta_bond); + state.truncate(); + } + } + Model::Heisenberg => { + for i in 0..n.saturating_sub(1) { + state.rxx(i, i + 1, theta_bond); + state.truncate(); + state.ryy(i, i + 1, theta_bond); + state.truncate(); + state.rzz(i, i + 1, theta_bond); + state.truncate(); + } + for i in 0..n { + state.rz(i, theta_site); + state.truncate(); + } + } + Model::Scramble => unreachable!("handled above"), + } + } +} + +/// `⟨0…0|O|0…0⟩` for TFIM; the `Z_0` coefficient for Heisenberg. +/// +/// The diagonal contraction is spelled out rather than delegated: `⟨0|Z|0⟩ = 1` +/// and `⟨0|X|0⟩ = ⟨0|Y|0⟩ = 0`, so exactly the X-free terms survive and each +/// contributes its own coefficient. +fn readout(state: &PauliSum>, n: usize, p: Params) -> f64 { + match p.model { + Model::Tfim => state + .iter() + .filter(|(word, _)| (0..n).all(|i| !word.get_xbit(i))) + .map(|(_, c)| *c) + .sum(), + Model::Heisenberg | Model::Scramble => state + .data() + .get(&PauliWord::from(site_word(n, 0, 'Z').as_str())) + .copied() + .unwrap_or(0.0), + } +} + +/// Run one model at width `n`, returning `(best seconds, final support, observable)`. +fn run(n: usize, p: Params) -> (f64, usize, f64) { + let base = seed::(n, p); + let gates = if p.model == Model::Scramble { + scramble_gates(n, p) + } else { + Vec::new() + }; + let mut best = f64::INFINITY; + let mut terms = 0usize; + let mut observable = f64::NAN; + for _ in 0..p.iters { + let mut state = base.clone(); + let t0 = Instant::now(); + propagate(&mut state, n, p, &gates); + best = best.min(t0.elapsed().as_secs_f64()); + terms = state.len(); + observable = readout(&state, n, p); + } + (best, terms, observable) +} + +/// Print the whole propagated support as `word coefficient`, largest first — +/// the format the driver diffs across engines. +fn dump(n: usize, p: Params) { + let mut state = seed::(n, p); + let gates = if p.model == Model::Scramble { + scramble_gates(n, p) + } else { + Vec::new() + }; + propagate(&mut state, n, p, &gates); + let mut out: Vec<(String, f64)> = state.iter().map(|(k, c)| (k.to_string(), *c)).collect(); + out.sort_by(|a, b| { + b.1.abs() + .partial_cmp(&a.1.abs()) + .unwrap() + .then(a.0.cmp(&b.0)) + }); + println!("# {} terms", out.len()); + for (word, coeff) in out { + println!("{word} {coeff:+.12e}"); + } +} + +macro_rules! dispatch { + ($f:ident, $n:expr, $p:expr) => { + match storage_bytes($n) { + 2 => $f::<2>($n, $p), + 4 => $f::<4>($n, $p), + 8 => $f::<8>($n, $p), + 16 => $f::<16>($n, $p), + 32 => $f::<32>($n, $p), + 64 => $f::<64>($n, $p), + b => panic!("no storage tier for {b} bytes (n = {})", $n), + } + }; +} + +fn main() { + let p = Params::from_env(); + let qubits: Vec = std::env::var("QUBITS") + .ok() + .map(|s| s.split(',').filter_map(|t| t.trim().parse().ok()).collect()) + .unwrap_or_else(|| vec![8, 12, 16, 20, 24, 28, 32]); + + if std::env::var("DUMP").is_ok() { + let n = qubits[0]; + dispatch!(dump, n, p); + return; + } + + eprintln!( + "ppvm {}: steps={} dt={} J={} h={} atol={:e} iters={}", + p.model.name(), + p.steps, + p.dt, + p.j, + p.h, + p.atol, + p.iters + ); + + println!("model,library,qubits,steps,dt,atol,time_s,terms,observable"); + for &n in &qubits { + let (t, terms, obs) = dispatch!(run, n, p); + println!( + "{},ppvm,{},{},{},{:e},{:.6},{},{:.12e}", + p.model.name(), + n, + p.steps, + p.dt, + p.atol, + t, + terms, + obs + ); + eprintln!(" n={n:3} {t:9.4}s {terms:>9} terms obs={obs:+.9e}"); + use std::io::Write; + std::io::stdout().flush().ok(); + } +} diff --git a/julia-benchmarks/Manifest.toml b/julia-benchmarks/Manifest.toml index 3c050de35..d846d81c7 100644 --- a/julia-benchmarks/Manifest.toml +++ b/julia-benchmarks/Manifest.toml @@ -1,8 +1,8 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.12.1" +julia_version = "1.12.6" manifest_format = "2.0" -project_hash = "2774cd547152d7fbb5c2095d2e68e179a1b41160" +project_hash = "52f0dd51f0e051696bfc24b6c265049f3b4a890d" [[deps.AbstractFFTs]] deps = ["LinearAlgebra"] @@ -208,6 +208,11 @@ git-tree-sha1 = "37ea44092930b1811e666c3bc38065d7d87fcc74" uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" version = "0.13.1" +[[deps.Combinatorics]] +git-tree-sha1 = "c761b00e7755700f9cdf5b02039939d1359330e1" +uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" +version = "1.1.0" + [[deps.Compat]] deps = ["TOML", "UUIDs"] git-tree-sha1 = "9d8a54ce4b17aa5bdce0ea5c34bc5e7c340d16ad" @@ -267,6 +272,12 @@ git-tree-sha1 = "9e2f36d3c96a820c678f2f1f1782582fcf685bae" uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" version = "1.9.1" +[[deps.Dictionaries]] +deps = ["Indexing", "Random", "Serialization"] +git-tree-sha1 = "a55766a9c8f66cf19ffcdbdb1444e249bb4ace33" +uuid = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4" +version = "0.4.6" + [[deps.Distances]] deps = ["LinearAlgebra", "Statistics", "StatsAPI"] git-tree-sha1 = "c7e3a542b999843086e2f29dac96a618c105be1d" @@ -307,7 +318,7 @@ version = "0.9.5" [[deps.Downloads]] deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" -version = "1.6.0" +version = "1.7.0" [[deps.EpollShim_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -478,6 +489,11 @@ git-tree-sha1 = "68c173f4f449de5b438ee67ed0c9c748dc31a2ec" uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a" version = "0.3.28" +[[deps.Indexing]] +git-tree-sha1 = "ce1566720fd6b19ff3411404d4b977acd4814f9f" +uuid = "313cdc1a-70c2-5d6a-ae34-0150d3930a38" +version = "1.1.1" + [[deps.Inflate]] git-tree-sha1 = "d1b1b796e47d94588b3757fe84fbf65a5ec4a80d" uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9" @@ -617,7 +633,7 @@ version = "0.6.4" [[deps.LibCURL_jll]] deps = ["Artifacts", "LibSSH2_jll", "Libdl", "OpenSSL_jll", "Zlib_jll", "nghttp2_jll"] uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" -version = "8.11.1+1" +version = "8.15.0+0" [[deps.LibGit2]] deps = ["LibGit2_jll", "NetworkOptions", "Printf", "SHA"] @@ -744,7 +760,7 @@ version = "1.11.0" [[deps.MozillaCACerts_jll]] uuid = "14a3606d-f60d-562e-9121-12d972cd8159" -version = "2025.5.20" +version = "2025.11.4" [[deps.MuladdMacro]] git-tree-sha1 = "cac9cc5499c25554cba55cd3c30543cff5ca4fab" @@ -812,7 +828,7 @@ version = "1.6.1" [[deps.OpenSSL_jll]] deps = ["Artifacts", "Libdl"] uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" -version = "3.5.1+0" +version = "3.5.4+0" [[deps.OpenSpecFun_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"] @@ -876,6 +892,20 @@ version = "0.7.3" [deps.PauliPropagation.weakdeps] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +[[deps.PauliStrings]] +deps = ["BitIntegers", "Combinatorics", "Dictionaries", "LinearAlgebra", "ProgressBars", "Random", "SparseArrays"] +git-tree-sha1 = "80ecc963d62d97541258fbc8fcdf2fc95a2e70f7" +uuid = "f07625cc-80e4-4099-a362-36a3484d9bcc" +version = "1.10.1" + + [deps.PauliStrings.extensions] + MathLinkPauliStringsExt = "MathLink" + SymbolicsPauliStringsExt = "Symbolics" + + [deps.PauliStrings.weakdeps] + MathLink = "18c93696-a329-5786-9845-8443133fa0b4" + Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" + [[deps.Pixman_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl"] git-tree-sha1 = "e4a6721aa89e62e5d4217c0b21bd714263779dda" @@ -885,7 +915,7 @@ version = "0.46.4+0" [[deps.Pkg]] deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "Random", "SHA", "TOML", "Tar", "UUIDs", "p7zip_jll"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -version = "1.12.0" +version = "1.12.1" weakdeps = ["REPL"] [deps.Pkg.extensions] @@ -963,6 +993,12 @@ git-tree-sha1 = "990016fb1508b0726a70039f39569720d054c78d" uuid = "efd6af41-a80b-495e-886c-e51b0c7d77a3" version = "0.1.7" +[[deps.ProgressBars]] +deps = ["Printf"] +git-tree-sha1 = "b437cdb0385ed38312d91d9c00c20f3798b30256" +uuid = "49802e3a-d2f1-5c88-81d8-b72133a6f568" +version = "1.5.1" + [[deps.PtrArrays]] git-tree-sha1 = "4fbbafbc6251b883f4d2705356f3641f3652a7fe" uuid = "43287f4e-b6f4-7ad1-bb20-aadabca52c3d" @@ -1580,9 +1616,9 @@ uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" version = "1.64.0+1" [[deps.p7zip_jll]] -deps = ["Artifacts", "Libdl"] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" -version = "17.5.0+2" +version = "17.7.0+0" [[deps.x264_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] diff --git a/julia-benchmarks/Project.toml b/julia-benchmarks/Project.toml index 2a0f4ed61..befd11060 100644 --- a/julia-benchmarks/Project.toml +++ b/julia-benchmarks/Project.toml @@ -2,4 +2,5 @@ BenchmarkPlots = "ab8c0f59-4072-4e0d-8f91-a91e1495eb26" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" PauliPropagation = "293282d5-3c99-4fb6-92d0-fd3280a19750" +PauliStrings = "f07625cc-80e4-4099-a362-36a3484d9bcc" StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" diff --git a/julia-benchmarks/benches/xbench_pp.jl b/julia-benchmarks/benches/xbench_pp.jl new file mode 100644 index 000000000..15087b28b --- /dev/null +++ b/julia-benchmarks/benches/xbench_pp.jl @@ -0,0 +1,129 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 +# +# The PauliPropagation.jl side of the cross-library Pauli-propagation +# benchmark. See `benchmarks/cross-library/README.md` for the shared circuit +# definitions, the parameter contract, and the CSV schema — every runner reads +# the same environment variables and prints the same columns. +# +# MODEL=tfim QUBITS=8,16,24 STEPS=10 DT=0.1 JCOUP=1.0 HFIELD=1.0 ATOL=1e-6 \ +# julia --project=@. -t1 benches/xbench_pp.jl + +using PauliPropagation +using Printf + +const MODEL = get(ENV, "MODEL", "tfim") +const STEPS = parse(Int, get(ENV, "STEPS", "10")) +const DT = parse(Float64, get(ENV, "DT", "0.1")) +const JCOUP = parse(Float64, get(ENV, "JCOUP", "1.0")) +const HFIELD = parse(Float64, get(ENV, "HFIELD", "1.0")) +const ATOL = parse(Float64, get(ENV, "ATOL", "1e-6")) +const ITERS = parse(Int, get(ENV, "ITERS", "3")) + +const THETA_BOND = 2 * JCOUP * DT +const THETA_SITE = 2 * HFIELD * DT + +"""The seed observable: `Σ_i Z_i` for TFIM, `Z_1` for Heisenberg.""" +function seed_operator(n::Int) + ps = PauliSum(n) + if MODEL == "tfim" + for i in 1:n + add!(ps, PauliString(n, [:Z], [i])) + end + else + add!(ps, PauliString(n, [:Z], [1])) + end + return ps +end + +"""One first-order Trotter step, gates in the order the shared spec fixes.""" +function trotter_step!(state, n::Int) + if MODEL == "tfim" + for i in 1:n + state = propagate(PauliRotation([:X], [i], THETA_SITE), state; min_abs_coeff = ATOL) + end + for i in 1:(n - 1) + state = propagate( + PauliRotation([:Z, :Z], [i, i + 1], THETA_BOND), state; min_abs_coeff = ATOL + ) + end + else + for i in 1:(n - 1) + for axes in ([:X, :X], [:Y, :Y], [:Z, :Z]) + state = propagate( + PauliRotation(axes, [i, i + 1], THETA_BOND), state; min_abs_coeff = ATOL + ) + end + end + for i in 1:n + state = propagate(PauliRotation([:Z], [i], THETA_SITE), state; min_abs_coeff = ATOL) + end + end + return state +end + +function run_model(n::Int) + state = seed_operator(n) + for _ in 1:STEPS + state = trotter_step!(state, n) + end + return state +end + +"""`⟨0…0|O|0…0⟩` for TFIM; the `Z_1` autocorrelator for Heisenberg.""" +function readout(state, n::Int) + if MODEL == "tfim" + return real(overlapwithzero(state)) + else + c = getcoeff(state, PauliString(n, [:Z], [1])) + return real(c) + end +end + +""" +Print the propagated support as `word coefficient`, largest first, site 1 +leftmost — the same format the Rust and Python runners emit under `DUMP=1` so +the driver can diff all four term-for-term. +""" +function dump_support(n::Int) + state = run_model(n) + terms = Dict{String,Float64}() + for (ps, c) in state + w = join([['I', 'X', 'Y', 'Z'][getpauli(ps, i) + 1] for i in 1:n]) + terms[w] = get(terms, w, 0.0) + real(c) + end + println("# $(length(terms)) terms") + for w in sort(collect(keys(terms)); by = w -> (-abs(terms[w]), w)) + println("$w $(Printf.@sprintf("%+.12e", terms[w]))") + end +end + +function main() + qubits = parse.(Int, split(get(ENV, "QUBITS", "8,12,16,20,24,28,32"), ",")) + if haskey(ENV, "DUMP") + dump_support(first(qubits)) + return + end + println( + stderr, + "PauliPropagation.jl $MODEL: steps=$STEPS dt=$DT J=$JCOUP h=$HFIELD atol=$ATOL " * + "iters=$ITERS threads=$(Threads.nthreads())", + ) + println("model,library,qubits,steps,dt,atol,time_s,terms,observable") + for n in qubits + # Warm up so the reported time excludes JIT for this width's type + # specialization. + st = run_model(n) + nterms = length(st) + obs = readout(st, n) + best = Inf + for _ in 1:ITERS + best = min(best, @elapsed run_model(n)) + end + println("$MODEL,pauli-propagation-jl,$n,$STEPS,$DT,$ATOL,$(round(best, sigdigits=7)),$nterms,$obs") + println(stderr, " n=$n $(round(best, digits=4))s $nterms terms obs=$obs") + flush(stdout) + end +end + +main() diff --git a/julia-benchmarks/benches/xbench_ps.jl b/julia-benchmarks/benches/xbench_ps.jl new file mode 100644 index 000000000..4a51d7380 --- /dev/null +++ b/julia-benchmarks/benches/xbench_ps.jl @@ -0,0 +1,163 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 +# +# The PauliStrings.jl side of the cross-library Pauli-propagation benchmark. +# See `benchmarks/cross-library/README.md` for the shared circuit definitions, +# the parameter contract, and the CSV schema. +# +# PauliStrings.jl's own front door for this is `evolve(H, O, tspan; +# method=Trotter())`, which builds the gate list from the Hamiltonian's internal +# string order. We build the `TrotterGate` list by hand instead, so the gate +# sequence is the one the shared spec fixes and the propagated operator is +# comparable term-for-term with the other three engines. +# +# MODEL=tfim QUBITS=8,16,24 STEPS=10 DT=0.1 JCOUP=1.0 HFIELD=1.0 ATOL=1e-6 \ +# julia --project=@. -t1 benches/xbench_ps.jl + +using PauliStrings +using Printf + +const MODEL = get(ENV, "MODEL", "tfim") +const STEPS = parse(Int, get(ENV, "STEPS", "10")) +const DT = parse(Float64, get(ENV, "DT", "0.1")) +const JCOUP = parse(Float64, get(ENV, "JCOUP", "1.0")) +const HFIELD = parse(Float64, get(ENV, "HFIELD", "1.0")) +const ATOL = parse(Float64, get(ENV, "ATOL", "1e-6")) +const ITERS = parse(Int, get(ENV, "ITERS", "3")) + +const THETA_BOND = 2 * JCOUP * DT +const THETA_SITE = 2 * HFIELD * DT + +"""The `PauliString` generator for a one- or two-site Pauli term.""" +function generator(n::Int, ops::Vector{String}, sites::Vector{Int}) + o = Operator(n) + if length(sites) == 1 + o += ops[1], sites[1] + else + o += ops[1], sites[1], ops[2], sites[2] + end + return o.strings[1] +end + +""" +The Trotter gate list in *application* order, per the shared spec. + +`trotter_step!` consumes gates in matrix-multiply order and applies them in +reverse, so the caller reverses this before handing it over. +""" +function gate_list(n::Int) + gates = PauliStrings.TrotterGate[] + if MODEL == "tfim" + for i in 1:n + push!(gates, PauliStrings.TrotterGate(generator(n, ["X"], [i]), THETA_SITE)) + end + for i in 1:(n - 1) + push!( + gates, + PauliStrings.TrotterGate(generator(n, ["Z", "Z"], [i, i + 1]), THETA_BOND), + ) + end + else + for i in 1:(n - 1) + for op in ("X", "Y", "Z") + push!( + gates, + PauliStrings.TrotterGate( + generator(n, [op, op], [i, i + 1]), THETA_BOND + ), + ) + end + end + for i in 1:n + push!(gates, PauliStrings.TrotterGate(generator(n, ["Z"], [i]), THETA_SITE)) + end + end + return gates +end + +"""The seed observable: `Σ_i Z_i` for TFIM, `Z_1` for Heisenberg.""" +function seed_operator(n::Int) + o = Operator(n) + if MODEL == "tfim" + for i in 1:n + o += "Z", i + end + else + o += "Z", 1 + end + return o +end + +function run_model(n::Int, gates) + o = seed_operator(n) + truncation(x) = cutoff(x, ATOL) + for _ in 1:STEPS + PauliStrings.trotter_step!(o, gates; truncation = truncation, truncate_every = 1) + end + return o +end + +"""`⟨0…0|O|0…0⟩` for TFIM; the `Z_1` autocorrelator for Heisenberg.""" +function readout(o::Operator, n::Int) + if MODEL == "tfim" + return real(expect(o, "0"^n)) + else + z1 = seed_operator(n) + return real(trace_product(z1, o; scale = 1)) + end +end + +""" +Print the propagated support as `word coefficient`, largest first, site 1 +leftmost — the same format the other three runners emit under `DUMP=1`. + +PauliStrings.jl carries `im^{#Y}` inside the stored coefficient (its `Matrix` +convention), so `op_to_strings` is the readout that puts the coefficients on the +same footing as the other engines' real ones. +""" +function dump_support(n::Int, gates) + o = run_model(n, gates) + coeffs, strings = op_to_strings(o) + terms = Dict{String,Float64}() + for (c, s) in zip(coeffs, strings) + w = replace(s, '1' => 'I') + terms[w] = get(terms, w, 0.0) + real(c) + end + println("# $(length(terms)) terms") + for w in sort(collect(keys(terms)); by = w -> (-abs(terms[w]), w)) + @printf("%s %+.12e\n", w, terms[w]) + end +end + +function main() + qubits = parse.(Int, split(get(ENV, "QUBITS", "8,12,16,20,24,28,32"), ",")) + if haskey(ENV, "DUMP") + n = first(qubits) + dump_support(n, reverse(gate_list(n))) + return + end + println( + stderr, + "PauliStrings.jl $MODEL: steps=$STEPS dt=$DT J=$JCOUP h=$HFIELD atol=$ATOL " * + "iters=$ITERS threads=$(Threads.nthreads())", + ) + println("model,library,qubits,steps,dt,atol,time_s,terms,observable") + for n in qubits + # `gate_list` is circuit construction, not propagation — built once and + # excluded from the timed region, as in every other runner. + gates = reverse(gate_list(n)) + # Warm up so the reported time excludes JIT. + o = run_model(n, gates) + nterms = length(o) + obs = readout(o, n) + best = Inf + for _ in 1:ITERS + best = min(best, @elapsed run_model(n, gates)) + end + println("$MODEL,pauli-strings-jl,$n,$STEPS,$DT,$ATOL,$(round(best, sigdigits=7)),$nterms,$obs") + println(stderr, " n=$n $(round(best, digits=4))s $nterms terms obs=$obs") + flush(stdout) + end +end + +main()