Metadata-Version: 2.4
Name: gbskernels
Version: 0.2.2
Summary: GPU-native, batched library of photonic sampling matrix functions (hafnian, loop hafnian, permanent, torontonian).
Author-email: Chaofan Zhou <chaofanzhou97@gmail.com>
License-Expression: Apache-2.0
Project-URL: Homepage, https://github.com/chaofanzhougit/GBSKernels
Project-URL: Documentation, https://github.com/chaofanzhougit/GBSKernels#readme
Project-URL: Issues, https://github.com/chaofanzhougit/GBSKernels/issues
Keywords: quantum,photonics,gaussian-boson-sampling,hafnian,permanent,torontonian,loop-hafnian,gpu,cuda,numerical-methods
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Scientific/Engineering :: Physics
Classifier: Topic :: Scientific/Engineering :: Mathematics
Requires-Python: >=3.12
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.26
Requires-Dist: mpmath>=1.3
Provides-Extra: publication
Requires-Dist: build==1.5.0; extra == "publication"
Requires-Dist: nanobind==2.13.0; extra == "publication"
Requires-Dist: numpy==2.4.6; extra == "publication"
Requires-Dist: mpmath==1.3.0; extra == "publication"
Requires-Dist: piquasso==8.0.1; extra == "publication"
Requires-Dist: python-flint==0.9.0; extra == "publication"
Requires-Dist: thewalrus==0.22.0; extra == "publication"
Dynamic: license-file

# GBSKernels

A GPU-native, batched library of the #P-hard matrix functions behind photonic
quantum sampling — the **permanent, hafnian, loop hafnian, and torontonian** —
with an explicit floating-point accuracy model that ranges from native double
precision to rigorous a-posteriori error bounds.

Gaussian boson sampling and its displaced- and threshold-detector variants, and
standard boson sampling, reduce computationally to evaluating *large batches of
independent, medium-sized* instances of these four functions. Each is an
alternating signed sum whose value can be many orders of magnitude smaller than
its terms, so it is cancellation-prone in floating point: both throughput and
*known* accuracy matter. GBSKernels targets exactly this workload — batched-first
kernels on CPU and CUDA, every function validated against independent
combinatorial ground truth, and a precision model in which the accuracy of every
result is either measured or provably bounded.

GBSKernels is designed to **complement**
[The Walrus](https://github.com/XanaduAI/thewalrus), the canonical CPU/C++
reference implementation, rather than to replace it: a GPU-native, batched
companion with an explicit accuracy characterization. The Walrus is used here as
one of several independent test oracles.

> **Status.** The CPU library is complete and `pip`-installable (pure Python,
> `numpy` + `mpmath` only). The v0.2.1 GitHub release includes a GPU validation
> manifest covering all 24 device gates, the Python binding smoke, adversarial
> and physical enclosure families, and the Jiuzhang Gate C probe. The session is
> release-commit-attributed through the environment-supplied `GBS_COMMIT`. The
> rsynced host had no Git metadata, so `tracked_dirty` is unknown; the evidence
> hashes and container digest do not hash-bind the uploaded source tree or loaded
> CUDA extension. Historical v0.2.0 RTX 4090 and A100 evidence is retained
> separately; see the [device evidence](results/README.md#v020-device-validation).
> This is pre-1.0 software and the public API may change.

Version 0.2.2 is a correctness follow-up: it charges the recursive-torontonian
double-double leaf accumulation against the operand magnitudes rather than the
post-addition result, so the returned enclosure radius stays sound under deep
single-step cancellation (returned values are unchanged). Version 0.2.1 added
subnormal-aware certified bounds and fail-closed certificate handling, made the
single-large real torontonian input contract explicit, removed inversion
roundoff skew from the xxpp construction, recorded whether artifacts came from
modified tracked sources, and packaged the corrected historical Jiuzhang
fixed-sample audit (exploratory, not a retrospective confirmatory result). See
the [changelog](CHANGELOG.md) for the release scope.

## Installation

The CPU package is available from
[PyPI](https://pypi.org/project/gbskernels/0.2.2/) and needs only `numpy` +
`mpmath`. The CUDA extension is a separate, optional source build (see
[`bindings/README.md`](bindings/README.md)):

```bash
python -m pip install gbskernels==0.2.2
```

## Quick start

```python
import numpy as np
import gbskernels

A = np.array([[1., 2.], [3., 4.]])
gbskernels.perm(A)                              # (10+0j)

# The API is batched-first: one call evaluates a stack of matrices.
rng = np.random.default_rng(0)
stack = rng.standard_normal((4096, 8, 8))
stack = stack + stack.transpose(0, 2, 1)        # symmetric, as the hafnian needs

gbskernels.haf_batched(stack)                   # double precision (default)
gbskernels.haf_batched(stack, precision="auto") # double precision + cancellation guard
gbskernels.haf_batched(stack, backend="gpu")    # CUDA, if the extension is built
gbskernels.haf(stack[0], precision="ref")       # arbitrary-precision reference
```

A short guided tour — closed-form checks, the `auto` tier catching real
double-precision cancellation, a GBS distribution, and a sampling run — runs
with wheel dependencies only:

```bash
python examples/gbs_demo.py
```

## Precision model

Precision is always an explicit tier, so the accuracy of a result is never left
implicit:

| Tier | Description |
|---|---|
| `"fp64"` (default) | Native double precision — the throughput path, with a measured accuracy boundary. |
| `"dd"` | Double-double arithmetic carried internally through the cancelling sum (a GPU tier); recovers a correct double-precision result where plain double precision cancels. |
| `"ref"` | Arbitrary-precision reference (`mpmath`); the accuracy ground truth. |
| `"auto"` | Double precision plus a per-evaluation cancellation indicator; evaluations flagged as risky are recomputed in a higher tier (`mpmath` on CPU, double-double on GPU). The indicator is a calibrated heuristic, not a certificate. |
| `"certified"` | The double-precision value together with a **rigorous a-posteriori error bound**: `\|value − exact\| ≤ abs_error_bound`, including explicit absolute terms for subnormal underflow (on the GPU, bound arithmetic uses per-instruction directed rounding). Non-finite inputs and invalid/non-finite certificates fail closed. Available for all four functions on CPU and GPU. With `rtol=`, a failed or refused certified-fp64 result escalates to the arbitrary-precision reference; GPU permanent and hafnian first try a certified double-double bound. The reference result is flagged as non-certified. `tor_single(..., dd=True)` provides a separate certified-double-double large-torontonian path. |

## Features

- **All four functions, batched, on CPU and CUDA.** Field-standard algorithms
  (Glynn/BB–FG for the permanent, the power-trace form for the hafnian and loop
  hafnian, subset-determinant inclusion–exclusion for the torontonian), each with
  an independent naive reference used for verification.
- **A measured accuracy boundary.** The double-precision–vs–double-double
  crossover is characterized for every function against an arbitrary-precision
  reference, on physical and adversarial inputs.
- **Certified evaluation.** Every value can be returned with a rigorous error
  bound whose enclosure of the true value is a hard test invariant — verified on
  closed forms, random ensembles, and adversarial cancellation families.
- **A device-resident workspace.** `gbskernels.Workspace` evaluates ragged
  batches with a persistent stream, pinned staging, and a zero-copy DLPack output
  that CuPy / PyTorch / JAX can consume without a device-to-host copy
  ([`docs/device_resident_contract.md`](docs/device_resident_contract.md)).
- **A conditional GBS sampler** (`sampling/`) that draws photon-number samples by
  the chain rule, evaluating each mode's batch of hafnians on the GPU, validated
  distributionally (total-variation distance below 0.03 against both the exact
  distribution and The Walrus, and a chi-square goodness-of-fit p > 10⁻³ against
  The Walrus).
- **Structure-aware kernels** for the sampling workload: a repeated-row
  finite-difference sieve for the loop hafnian, and a recursive prefix-Cholesky
  torontonian — including a single-large mode that splits one evaluation across
  the GPU to dimension 64 (32 modes).
- **A fail-closed confirmatory workflow** for the public Jiuzhang example,
  covering exposure audit, frozen design, future-beacon selection,
  content-addressed evaluation, refusal recovery, simultaneous coherence-grid
  inference, predictive checks, and release verification
  ([protocol](docs/confirmatory_v2.md),
  [tools](examples/jiuzhang/README.md)). The repository ships templates and
  validation contracts, not a fabricated registration or completed v2 outcome.

## Verification

Verification is treated as the deliverable ([`docs/DESIGN.md`](docs/DESIGN.md), §8).
Five layers are exercised by the test suite on every commit:

1. **Independent combinatorial ground truth** — perfect-matching counts and
   closed forms (`perm(J_n) = n!`, `haf(J_{2n}) = (2n−1)!!`, …), sharing no code
   with any existing library.
2. **Differential testing against The Walrus** across sizes, seeds, and
   real/complex inputs.
3. **Property-based invariants** (Hypothesis): scaling, permutation, transpose,
   the loop-hafnian → hafnian reduction, and others.
4. **Statistical / end-to-end validation** — the kernels compute GBS
   distributions that sum to one and match The Walrus per pattern in all three
   detector regimes.
5. **Numerical-accuracy characterization** against the arbitrary-precision
   reference, including adversarial cancellation families that defeat double
   precision and must survive double-double.

The Jiuzhang v2 tooling adds workflow-contract tests for canonical hashing,
registration timing, exclusion-ledger completeness, selection, immutable run
reduction, refusal recovery, inference, reconstruction, and release assembly.

To run the public v0.2.1 GPU validation from a fresh clone, first stage the
small hash-bound inputs (the raw USTC and Zenodo archives remain
external; see [`THIRD_PARTY_NOTICES.md`](THIRD_PARTY_NOTICES.md) for data
attribution and licensing scope):

```bash
python scripts/prepare_validation_data.py
```

The development discipline is CPU-first: host-shim-compatible CUDA kernels are
compiled and run through the CPU pre-flight (part of the test suite), then the
full source set is checked on-device against independent CPU references before
any throughput result is considered publishable. Device-only variants such as
the warp-specialized permanent are covered by the on-device gates; profiler
diagnostics collected before those gates remain provisional until they pass.

## Publication evidence workflow (prospective)

The repository also provides a stricter publication-session harness. It is
separate from the v0.2.1 evidence described above: no GPU result from this new
harness is being claimed here, and its artifacts become evidence only after an
actual run completes all fail-closed checks.

The on-box driver consumes a release archive and a separate immutable extraction,
checks both SHA-256 identities before and after the run, and requires the full Git
commit, Git tree, and digest-pinned container identity. Builds, virtual
environment, wheel, CUDA extension, logs, and evidence are created outside the
source tree. The resulting manifest binds the release source to compiler
commands, tool versions, the embedded-cubin inventory and SASS dump, available
PTX, the loaded extension, the built wheel, device
metadata, dependency resolution, and every retained artifact.

After the 24 device gates, the registered science profile is configured to run:

- an independent python-flint/Arb enclosure campaign with 312 dense small cases
  and one structured quarter-identity case at each `k = 25,...,32`; the latter
  probes the production recursion and certificate at those depths, but is not
  general dense large-`k` validation; and
- a frozen same-input torontonian benchmark through 20 modes, comparing the
  GBSKernels double-double GPU path with The Walrus and Piquasso recursive CPU
  implementations on 15 physical real-quadrature lossy-state matrices. The
  prospective schema-v4 artifact gives every exact frozen binary64 matrix an
  independent dense python-flint/Arb interval, evaluated outside the timed
  region, and checks whether the GBSKernels-reported DD radius contains that
  interval. The Walrus and Piquasso fp64 values carry no claimed error bounds.
  Pairwise tolerance flags are descriptive only: they neither accept nor reject
  the artifact nor filter timing rows. Because the timings compare GPU DD with
  CPU fp64 implementations, they do not establish a causal hardware or
  algorithmic speedup.

The Vast.ai launcher adds a generated archive adapter, digest-pinned image,
deterministic offer validation, dry-run and explicit-spend modes, cost/lifetime
ceilings, bounded ambiguous-create recovery and SSH handling, create-only
retrieval, receipts, and exact-ID non-interactive teardown with provider-list
absence verification on normal, failure, and handled-signal paths. See the
[`scripts/` workflow](scripts/README.md),
[`bench/` protocol](bench/README.md),
[`envs/` dependencies](envs/README.md), and
[`Jiuzhang Arb campaign`](examples/jiuzhang/README.md#publication-arb-enclosure-campaign-prospective)
for the precise scope and caveats.

## Performance

The library is batched-first: throughput comes from evaluating thousands of
independent instances per launch, and the benchmark harness reports where the
CPU wins (small single evaluations) as plainly as where the GPU wins (large
batches), always alongside the accuracy each point was measured at. The
structure-aware kernels give substantial speedups in their regimes — the
recursive torontonian over the LU form, and the finite-difference sieve over the
expanded power-trace at repeated-row patterns. Raw benchmark data is written,
append-only, to `results/`; the protocol is described in
[`docs/benchmark_protocol.md`](docs/benchmark_protocol.md).

```bash
python -m bench.accuracy                          # the accuracy boundary, all four functions
python -m bench.throughput --func perm --sizes 4,6,8,10
bash scripts/gpu_session.sh                       # full on-device session (on a CUDA host)
```

## Scope and limitations

These are correctness-first kernels with a deliberately narrow scope:

- **One evaluation per GPU thread** is the default mapping. Warp/block-cooperative
  variants exist and are dispatched only where they measurably win (the
  cooperative permanent); for the hafnian, loop hafnian, and torontonian the
  per-thread linear algebra is memory-bound, and the cooperative strategy does
  not help them.
- **The four functions do not share a single subset-enumeration engine.** Only
  the permanent uses the Gray-code delta walk; the others enumerate masks
  independently, because their per-subset work is memory-bound.
- **Fixed size limits** (per-thread buffers): permanent ≤ 28, hafnian ≤ 20 (DD
  ≤ 16), loop hafnian ≤ 20 (DD ≤ 14), torontonian 2n ≤ 24. Larger inputs raise on
  `backend="gpu"`; use the CPU backend. The single-large recursive torontonian
  extends to dimension 64.
- **The double-double tier is internal**: kernels carry double-double through the
  cancelling summation and collapse to `complex128` on output — they recover a
  correct double-precision answer under cancellation rather than exposing
  extended-precision results. The double-double torontonian is real-input only.
- **The sampler is pure-state and hybrid host-orchestrated** (the host drives the
  chain rule; the GPU evaluates each step's batch of hafnians). A fully
  on-device variant is available but limited to double precision within a
  submatrix-size cap.
- **Single GPU, CUDA only** — no multi-GPU, ROCm, or automatic differentiation.

## Repository layout

```
gbskernels/     batched-first Python API (the public surface)
cpu_ref/        independent CPU reference implementations (double + double-double)
highprec_ref/   arbitrary-precision references (mpmath; publication-only Arb)
core/           CUDA C++: subset-enumeration utilities, the kernels, on-device gates
bindings/       nanobind extension exposing the GPU backend to Python
sampling/       boson-sampling orchestration and the conditional GBS sampler
examples/       runnable CPU demo plus Jiuzhang reproduction and v2 workflow tools
tests/          the five-layer verification suite
bench/          accuracy and throughput harnesses
scripts/ envs/  scripted GPU-session runner and container definitions
docs/           design, benchmark/device contracts, and v2 protocol/templates
results/        append-only validation, gate, benchmark, and profiler artifacts
```

## Development

```bash
uv sync                      # Python 3.12 environment, dependencies, editable install
uv run pytest                # the CPU verification suite
uv run pytest -m layer1      # the independent combinatorial ground truth
uv run pytest -m "not slow"  # the fast tier
```

Building the CUDA extension requires `nvcc` and is a separate, optional step
([`bindings/README.md`](bindings/README.md)); the pure-Python wheel never
requires it. The full design rationale — algorithms, precision strategy, and the
verification contract — is in [`docs/DESIGN.md`](docs/DESIGN.md).

## Relationship to The Walrus

GBSKernels would not exist without [The Walrus](https://github.com/XanaduAI/thewalrus),
which remains the canonical reference implementation and this library's primary
differential oracle. The scope here is deliberately narrower: the four batched
kernels and their precision characterization.

## License

[Apache-2.0](LICENSE).
