diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
new file mode 100644
index 0000000..711302a
--- /dev/null
+++ b/.github/workflows/ci.yml
@@ -0,0 +1,122 @@
+name: CI
+
+on:
+ push:
+ branches: [main]
+ pull_request:
+
+jobs:
+ check:
+ runs-on: ubuntu-latest
+ strategy:
+ fail-fast: false
+ matrix:
+ python-version: ["3.10", "3.11", "3.12"]
+ steps:
+ - uses: actions/checkout@v4
+
+ - uses: actions/setup-python@v5
+ with:
+ python-version: ${{ matrix.python-version }}
+ cache: pip
+
+ - name: Install
+ run: |
+ python -m pip install --upgrade pip
+ pip install torch --index-url https://download.pytorch.org/whl/cpu
+ pip install -e ".[dev]"
+
+ - name: Lint
+ run: ruff check src tests scripts
+
+ - name: Type check
+ run: mypy
+
+ - name: Test
+ run: pytest tests/ -q
+
+ # The playable page quotes measured figures. This fails if any of them
+ # has drifted from what the artefacts in the repository actually say.
+ - name: Fail if the numbers on the playable page are stale
+ if: matrix.python-version == '3.12'
+ run: python scripts/site_facts.py --check
+
+ browser-physics:
+ # The browser runs a hand port of src/cueai/physics/. This regenerates the
+ # reference outcomes from the Python simulator and replays them in Node, so
+ # the two implementations cannot drift apart without the build going red.
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v4
+
+ - uses: actions/setup-python@v5
+ with:
+ python-version: "3.12"
+ cache: pip
+
+ - uses: actions/setup-node@v4
+ with:
+ node-version: "20"
+
+ - name: Install
+ run: |
+ python -m pip install --upgrade pip
+ pip install -e .
+
+ - name: Fail if the committed reference shots no longer match the simulator
+ run: python scripts/export_parity_cases.py --check
+
+ - name: Check the browser physics against the reference
+ run: node web/test/parity.mjs --verbose
+
+ - name: Play the bot against itself
+ run: node web/test/selfplay.mjs --games 6 --difficulty club
+
+ - name: Load the page in Chrome and play a game
+ run: |
+ npm install --no-save puppeteer-core
+ python -m http.server 8123 --directory web &
+ npx --yes wait-on http://localhost:8123/index.html
+ node web/test/browser.mjs --url http://localhost:8123/index.html --games 1
+
+ # The layer above: pointer capture, the drag threshold, and which element
+ # owns the spacebar. Driven with a real cursor, because calling the
+ # handlers directly is exactly what would miss a regression here.
+ - name: Drive the page with a real cursor and keyboard
+ run: node web/test/input.mjs --url http://localhost:8123/index.html
+
+ pipeline:
+ # Proves the data -> train -> export -> serve path still works end to end,
+ # on a dataset small enough to finish in a couple of minutes.
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v4
+
+ - uses: actions/setup-python@v5
+ with:
+ python-version: "3.12"
+ cache: pip
+
+ - name: Install
+ run: |
+ python -m pip install --upgrade pip
+ pip install torch --index-url https://download.pytorch.org/whl/cpu
+ pip install -e ".[dev]"
+
+ - name: Train on a small dataset
+ run: python -m cueai.ml.train --n-samples 400 --epochs 10
+
+ - name: Benchmark
+ run: python scripts/benchmark.py --repeats 50
+
+ - name: Render figures
+ run: python scripts/make_figures.py
+
+ - uses: actions/upload-artifact@v4
+ with:
+ name: cueai-pipeline-output
+ path: |
+ models/metrics.json
+ models/latency.json
+ docs/BENCHMARKS.md
+ docs/assets/*.png
diff --git a/.github/workflows/pages.yml b/.github/workflows/pages.yml
new file mode 100644
index 0000000..fdc57d8
--- /dev/null
+++ b/.github/workflows/pages.yml
@@ -0,0 +1,66 @@
+name: Deploy the playable table
+
+# The game is dependency-free ES modules, so "building" it is copying web/.
+# It still has to pass the physics parity check before it goes live: a page
+# that plays differently from the reference simulator would undercut the point
+# of publishing it.
+
+on:
+ push:
+ branches: [main]
+ workflow_dispatch:
+
+permissions:
+ contents: read
+ pages: write
+ id-token: write
+
+concurrency:
+ group: pages
+ cancel-in-progress: true
+
+jobs:
+ verify:
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v4
+
+ - uses: actions/setup-python@v5
+ with:
+ python-version: "3.12"
+ cache: pip
+
+ - uses: actions/setup-node@v4
+ with:
+ node-version: "20"
+
+ - name: Install the reference simulator
+ run: |
+ python -m pip install --upgrade pip
+ pip install -e .
+
+ - name: Check the browser physics against the reference
+ run: |
+ python scripts/export_parity_cases.py --check
+ node web/test/parity.mjs
+
+ - name: Play the bot against itself
+ run: node web/test/selfplay.mjs --games 4 --difficulty club
+
+ deploy:
+ needs: verify
+ runs-on: ubuntu-latest
+ environment:
+ name: github-pages
+ url: ${{ steps.deployment.outputs.page_url }}
+ steps:
+ - uses: actions/checkout@v4
+
+ - uses: actions/configure-pages@v5
+
+ - uses: actions/upload-pages-artifact@v3
+ with:
+ path: web
+
+ - id: deployment
+ uses: actions/deploy-pages@v4
diff --git a/.gitignore b/.gitignore
index b79eb94..02134a7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -15,3 +15,8 @@ data/processed/
*.log
.DS_Store
.cppbuild/
+
+# Puppeteer, installed on demand by web/test/browser.mjs
+node_modules
+package-lock.json
+docs/assets/.frames/
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..8bbc8b5
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2026 CueAI contributors
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..d4e8240
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,77 @@
+.PHONY: help setup test lint typecheck check train bench figures facts api ui play \
+ parity parity-check selfplay browser input capture web clean all
+.DEFAULT_GOAL := help
+
+PY ?= python3
+NODE ?= node
+# Defaults reproduce the published numbers in models/metrics.json.
+SAMPLES ?= 20000
+EPOCHS ?= 300
+GAMES ?= 20
+PORT ?= 8123
+
+help: ## Show available targets
+ @grep -E '^[a-z-]+:.*?## ' $(MAKEFILE_LIST) | awk -F':.*## ' '{printf " \033[36m%-10s\033[0m %s\n", $$1, $$2}'
+
+setup: ## Install the package with dev extras
+ $(PY) -m pip install -e ".[dev]"
+
+test: ## Run the test suite
+ $(PY) -m pytest tests/ -q
+
+lint: ## Check style
+ $(PY) -m ruff check src tests scripts
+
+typecheck: ## Run static type checks
+ $(PY) -m mypy
+
+check: lint typecheck test ## Everything CI runs on the Python side
+
+parity: ## Re-export reference shots and check the browser physics against them
+ $(PY) scripts/export_parity_cases.py
+ $(NODE) web/test/parity.mjs --verbose
+ $(PY) scripts/site_facts.py
+
+parity-check: ## Verify the committed reference shots without rewriting them
+ $(PY) scripts/export_parity_cases.py --check
+ $(NODE) web/test/parity.mjs
+
+selfplay: ## Play the bot against itself headlessly, exercising every rule
+ $(NODE) web/test/selfplay.mjs --games $(GAMES) --difficulty sharp --vs relaxed
+
+browser: ## Load the page in Chrome and play a game (needs puppeteer-core)
+ $(NODE) web/test/browser.mjs --url http://localhost:$(PORT)/index.html --games 2
+ $(NODE) web/test/input.mjs --url http://localhost:$(PORT)/index.html
+
+capture: ## Re-record the screenshots and the clip in the README
+ $(NODE) web/test/capture.mjs --url http://localhost:$(PORT)/index.html
+
+web: parity selfplay ## Every check that does not need a browser
+
+play: ## Serve the game at http://localhost:$(PORT)
+ @echo "CueAI is at http://localhost:$(PORT)"
+ @cd web && $(PY) -m http.server $(PORT)
+
+train: ## Generate data and train the residual model
+ $(PY) -m cueai.ml.train --n-samples $(SAMPLES) --epochs $(EPOCHS)
+
+bench: ## Measure latency and rewrite docs/BENCHMARKS.md
+ $(PY) scripts/benchmark.py
+
+figures: ## Render the README figures into docs/assets
+ $(PY) scripts/make_figures.py
+
+facts: ## Rewrite the numbers the playable page quotes
+ $(PY) scripts/site_facts.py
+
+api: ## Serve the prediction API on :8000
+ $(PY) -m uvicorn cueai.api.main:app --reload --port 8000
+
+ui: ## Launch the desktop table (needs the ui extra)
+ $(PY) -m cueai.ui.app
+
+all: train bench figures facts ## Reproduce every published number and figure
+
+clean: ## Remove generated artefacts
+ rm -rf .pytest_cache .mypy_cache .ruff_cache **/__pycache__
+ rm -f models/*.pt models/*.onnx models/*.joblib
diff --git a/README.md b/README.md
index 8344d61..5cf97bf 100644
--- a/README.md
+++ b/README.md
@@ -1,57 +1,493 @@
# CueAI
-Physics-informed AI simulation for billiards that combines classical mechanics with
-regression models to predict realistic ball trajectories from spin, launch angle,
-velocity, cushion interactions, and table surface variations.
+[](https://github.com/BruceMoseti/cueai/actions/workflows/ci.yml)
+
+
-| Layer | Tech |
-|-------|------|
-| Physics core | C++ (optional) + Python NumPy simulator |
-| ML | PyTorch → ONNX, Scikit-learn baselines |
-| Data | Pandas, NumPy synthetic shot datasets |
-| Vision | OpenCV table / cue overlay helpers |
-| Backend | FastAPI |
-| Frontend | PyQt6 interactive table |
+**A physics simulator for billiards, a closed-form solution that replaces it, and
+a learned model that corrects what the closed form misses — roughly 7,600x
+faster than integration, with the accuracy measured and the failure mode stated.
+Plus a browser game against a search-based bot, running the same physics.**
-## Quick start
+Predicting where the balls come to rest costs **about 4.6 seconds per rack shot**
+by numerical integration. This project reduces that to **0.60 ms** with a mean
+error of **98 mm for direct shots** on a 2.54 x 1.27 m table, and — the part that
+makes it usable — tells you in advance which of its own predictions to trust.
+
+Every number below comes from this repository: `make all` regenerates the tables
+and figures from scratch, and `make test` checks the physics claims against their
+closed-form references.
+
+---
+
+## Play it
+
+### **[▶ Play eight-ball against the bot](https://brucemoseti.github.io/cueai/)**
+
+
+
+*Recorded by `web/test/capture.mjs` driving the real page in Chrome, at the
+"Quick" playback the page offers, with the bot's search pauses capped at
+0.8 s. The shots are the game's, not an animation.*
+
+No install, no build step, no framework — `web/` is plain ES modules, and the
+whole thing is served as static files. Aim with the mouse, drag back from the
+cue ball to strike, click the cue-ball diagram to move the tip off centre. Every
+shot is judged by a full eight-ball rules engine, and the log under the table
+keeps its reasoning on screen: which ball was struck first, how many rails were
+found afterwards, and which of those facts made a shot a foul.
+
+
+
+The browser is not running a lookalike physics engine. `web/js/physics.js` is a
+hand port of `src/cueai/physics/`, and the port is measured rather than
+asserted: `scripts/export_parity_cases.py` runs 35 shots through the Python
+simulator — draw, follow, english off two rails, thin cuts, clusters and full
+sixteen-ball breaks — and `web/test/parity.mjs` replays every one of them in
+Node and compares where each ball stopped.
+
+| | |
+|---|---:|
+| Reference shots replayed | 35 |
+| Worst disagreement, any ball, any shot | **1.2 × 10⁻³ mm** |
+| Table time compared | 156 s |
+| Browser against the Python reference | **65× faster** |
+
+Continuous integration runs that check on every push and refuses to deploy the
+page if the two drift apart.
+
+Twenty headless bot-against-bot games walk every branch of the rules, and check
+on each of the roughly eight hundred thousand physics steps they take that no
+two balls are ever sharing space and nothing is ever inside a cushion. The worst
+overlap seen is 0.5 mm on a 57 mm ball — a fifth of a screen pixel — which is
+the difference between "the solver converges" and "it looks like it does". They
+also assert that each break moves at least ten of the fifteen; it moves twelve,
+and the check exists because a break that clips the apex leaves the triangle
+standing while looking, convincingly, like a limitation of the physics.
+
+**The opponent is a search, and that is the argument.**
+
+
+
+Aiming needs no learned model: the ghost-ball construction is exact, and a test
+asserts that aiming at it pots the ball while half a degree either side misses.
+So the bot spends nothing on aiming. It enumerates every ball-and-pocket pair in
+closed form, discards what is blocked or cut too thin, and spends its entire
+budget *simulating* the survivors to see what each one leaves behind.
+
+It then reports what that cost, every turn, in its own units: four pot lines
+solved in closed form, sixteen futures simulated, a minute of table time inside
+a third of a second of yours. That is what cheap physics buys, and it is why the
+strength setting is a rollout budget rather than an adjective.
+
+
+
+**The panels are there to be read.**
+
+
+
+A live inspector traces the cue ball's centre-of-mass speed
+against its contact-point slip while the shot is in flight, and draws the
+predicted `5/7·v₀` rolling speed as a line the measured curve has to land on.
+The yellow slip curve collapsing to zero exactly where the green speed curve
+flattens onto the blue line is the cloth model's central claim, happening in
+front of you, at whatever speed you set the playback to. It is a prediction from
+the mechanics, not a fitted parameter.
+
+The prediction is drawn only across the stretch of the shot it is a prediction
+about. Hit a rail after the handover, as in the trace above, and the line stops
+there, because the speed on the far side of a collision is that collision's
+business. Hit something *before* the handover and the panel withdraws the line
+entirely and says why, rather than showing the simulation missing a target it
+was never aiming at. The time axis expands the first fraction of a second when
+the tail is long, because the handover is over in about 150 ms and everything
+else is a ball rolling.
+
+
+
+**And there is prose under the table.**
+[The explainer](https://brucemoseti.github.io/cueai/#how) covers the cloth
+model, the parity harness, the bot's search, where the learned surrogate helps
+and where it does not, and
+[the multi-ball contact bug](#the-bug-the-tests-could-not-see) the single-ball
+validation suite could never have caught. Every figure it quotes is written into
+it by `scripts/site_facts.py` from the artefacts of an actual run, and CI fails
+if the page and the measurements disagree.
```bash
-cd ~/Projects/cueai
-python3 -m venv .venv
-source .venv/bin/activate
-pip install -r requirements.txt
+make play # serve it at http://localhost:8123
+make web # parity against Python, then 20 headless bot-vs-bot games
+make browser # play two games in Chrome, then 15 real-cursor interaction checks
+```
+
+That last one is worth a sentence, because it is the layer most browser tests
+skip. Calling the page's own functions proves the modules wire together; it says
+nothing about pointer capture, the drag threshold that separates lining a shot
+up from playing it, or which element owns the spacebar. `web/test/input.mjs`
+moves an actual cursor and presses actual keys, and asserts the promises the
+interface makes to a player: a quick click does not shoot, drawing the cue back
+sets the power without disturbing the aim, shift eases the aim to within a
+hundredth of a degree, the spin widget clamps at the miscue limit, and a
+dropdown that has focus keeps its own spacebar.
+
+---
+
+## Results
+
+Held-out test set of 4,000 shots from 20,000 simulated shots. Error is the
+distance between the predicted and the simulated resting position, averaged over
+the cue ball and the object ball. Full tables in
+[docs/BENCHMARKS.md](docs/BENCHMARKS.md).
+
+| Method | Cost per shot | Mean error | Direct shots (no cushion) | R² |
+|---|---:|---:|---:|---:|
+| Numerical simulator, 16 balls | 4.6 s | — (ground truth) | — | — |
+| Closed-form solver, no fitting | 0.27 ms | 494 mm | 114 mm | 0.01 |
+| Gradient boosting on the same features | 2.6 ms | 382 mm | 162 mm | 0.52 |
+| **Closed form + learned residual** | **0.60 ms** | **376 mm** | **98 mm** | 0.43 |
+
+The learned residual has the lowest error of the three, is four times cheaper to
+evaluate than the boosted trees, and is the only one that improves on the physics
+for the shots where physics is nearly sufficient. Gradient boosting posts the
+higher R² by hedging toward the middle of the table on shots nobody can predict,
+which flatters the variance-explained metric and costs it 64 mm on the shots that
+matter.
+
+### What that average hides
+
+Two thirds of sampled shots never actually reach the object ball. It stays where
+it started, which the closed-form baseline predicts *exactly*, so those rows
+donate a free zero to half of the error metric. Splitting them out, in mm:
+
+| Ball-ball contact | Share | Closed form cue / object | Boosting cue / object | CueNet cue / object |
+|---|---:|---:|---:|---:|
+| no | 67.6% | 585 / **0** | 457 / 77 | 411 / 16 |
+| yes | 32.4% | 1033 / 798 | **638 / 610** | 736 / 696 |
+
+Two things fall out of this that are worth saying plainly.
+
+The residual formulation earns its keep on the shots where nothing happens: it
+adds 16 mm of spurious object-ball motion against gradient boosting's 77 mm,
+because "predict zero correction" is its default rather than something it has to
+learn.
+
+And it loses on the shots where a collision has to be modelled: 736 mm against
+638 mm for the cue ball. That is the cost of anchoring to a baseline with a blind
+spot — the closed-form solver has no ball-ball contact model at all, so on a third
+of the shot space the residual is asked to undo a metre of error rather than
+refine a good guess. Giving the closed form even a crude ghost-ball collision
+would be the highest-value next change, and it is a physics change, not an ML one.
+
+### Where it stops working, and why that is the interesting part
+
+
+
+A resting position is a smooth function of the shot until the ball starts
+ricocheting between cushions. After two or three rail contacts, a millimetre of
+cue placement moves the outcome by a table length. The error breakdown above is
+published instead of hidden inside a single average, because "this model predicts
+direct and one-rail shots to about 10 cm and multi-rail scatter not at all" is a
+usable statement, while "376 mm mean error" is not.
+
+### Knowing which predictions to trust, before making them
+
+
+
+The breakdown above is sliced by what the simulator *did*, which you only know
+after paying the 4.6 seconds. That makes it an autopsy, not a control.
+
+But the closed-form solver reports how many cushions it *expects* on the way to
+its answer, and that number is already computed as part of the prediction, so it
+is free. It turns out to be a good enough proxy for "is this shot chaotic" to use
+as a gate:
+
+| Answer only when the solver expects | Coverage | Mean error |
+|---|---:|---:|
+| No cushion | 9.8% | **100 mm** |
+| At most one cushion | 27.5% | 189 mm |
+| At most two cushions | 50.0% | 251 mm |
+| Anything (no gate) | 100% | 376 mm |
+
+So the fast path is not a 376 mm model. It is a 100 mm model that knows it should
+decline nine shots in ten, or a 251 mm model over half the shot space, and the
+shots it declines can be sent to the simulator. A surrogate that reports its own
+applicability domain can be deployed; one with a single headline error cannot.
+The learned residual is also the most accurate of the three models at every
+coverage level, which is a stronger claim than its 1.6% edge on the overall mean.
+
+### Speed
+
+
+
+## The physics is verified, not asserted
-# Generate data + train ML model
-python -m cueai.ml.train --n-samples 4000 --epochs 40
+
-# API server
-uvicorn cueai.api.main:app --reload --port 8000
+The simulator is tested against closed-form solutions and conservation laws, not
+against its own earlier output. A struck ball must begin rolling at exactly
+`5/7` of its launch speed after sliding `12v₀²/(49 μ_s g)`; frictional collisions
+must conserve momentum to machine precision; no contact may create energy. See
+[docs/VALIDATION.md](docs/VALIDATION.md) for the full table of properties,
+tolerances and measured deviations — and for the list of effects deliberately not
+modelled, such as cue ball squirt and swerve.
-# Desktop UI
-python -m cueai.ui.app
+This mattered. The original implementation had the contact-point slip velocity
+and the friction torque in opposite handedness, so friction drove a struck ball
+*away* from rolling: it slid until it stopped, the rolling speed was `0` instead
+of `5/7 v₀`, and stopping distances were four times short. Nothing in the code
+looked wrong. A closed-form comparison found it immediately.
+
+
+
+*One stroke speed, three cue tip heights. Backspin brings the cue ball back
+behind where it started, a centre-ball hit stops it dead at the object ball,
+topspin sends it through. All three come out of the cloth model, not from
+special-casing.*
+
+### The bug the tests could not see
+
+Closed-form checks exercise one ball at a time, so a defect that only exists
+*between* balls survives all of them. This one did.
+
+The collision solver counted two balls as touching when the gap between their
+surfaces was under `1e-4 m`, and `rack.py` built the triangle with a `1e-4 m`
+clearance. All thirty contacts in the rack therefore sat exactly on the
+threshold, and which side each one fell on came down to whether `hypot` rounded
+up or down. Sixteen registered. Fourteen did not.
+
+The consequence was a break propagating through a contact graph with holes in
+it: balls in the middle of the rack came out of a full-power break having barely
+moved, and — the tell — the table opened up *less* the harder it was struck. No
+test failed. It was found by measuring break spread against cue speed and
+getting the sign wrong.
+
+The fix was to give the tolerance a name (`CONTACT_BAND = 1e-5`), use that one
+everywhere the solver asks whether two balls are in contact, and rack the balls
+actually touching so nothing sits on the boundary. Three regression tests now
+hold it: every contact in a fresh rack is inside the band, the band is orders of
+magnitude away from both the floating-point noise below it and the rack spacing
+above it, and resolving an untouched rack changes nothing and converges on the
+first pass.
+
+The property that caught it is now a fourth test, because it is exactly the kind
+nobody writes down: mean distance from the centre of the pack has to rise
+strictly from 3 to 6 to 9 m/s. At 8.2 m/s, mean pair separation after a break
+went from 0.478 m before the fix to 0.637 m after, and it now increases with
+every increase in cue speed instead of falling.
+
+The general point is that a validation suite is only as broad as the situations
+it constructs. Ten exact single-ball tests passing is not evidence about sixteen
+balls in contact.
+
+
+
+## Try it
+
+```bash
+git clone https://github.com/BruceMoseti/cueai && cd cueai
+python -m venv .venv && source .venv/bin/activate
+pip install -e ".[dev]"
+
+make test # the whole suite, including closed-form validation
+make check # ruff, mypy, tests — what CI runs on the Python side
+make play # serve the playable table at http://localhost:8123
+make web # browser physics against Python, then 20 headless games
+make all # regenerate the dataset, model, benchmarks and figures
+```
+
+`make all` takes about 20 minutes end to end on eight cores, of which 15 are
+simulating the 20,000 training shots. Generation is seeded per sample, so the
+dataset it produces is identical to the one behind the numbers above.
+
+### Serve predictions
+
+```bash
+make api # http://localhost:8000/docs
```
-## Architecture
+```bash
+# Sub-millisecond estimate: closed form plus learned residual
+curl -s localhost:8000/predict/fast -H 'content-type: application/json' \
+ -d '{"speed": 2.2, "angle_deg": 8, "english_y": -0.3, "cue_x": 0.6, "cue_y": 0.6}'
+# Full simulation, with the fast estimate alongside and the gap between them
+curl -s localhost:8000/predict -H 'content-type: application/json' \
+ -d '{"speed": 2.2, "angle_deg": 8, "full_rack": false, "obj_x": 1.4, "obj_y": 0.7}'
```
-Shot params (V, θ, spin ω, table μ)
- │
- ▼
-┌───────────────────┐ residual correction
-│ Physics simulator │ ──► ┌─────────────────┐
-│ (cloth, collide, │ │ CueNet (PyTorch)│──► ONNX
-│ cushions) │ └─────────────────┘
-└───────────────────┘ │
- │ ▼
- └──────────► fused trajectory ──► API / PyQt UI
+
+### Interactive table
+
+```bash
+pip install -e ".[ui]"
+make ui
```
-## Resume bullets (use the AI-focused version)
+Drag balls, aim, shoot. Each rack shot runs the reference simulator, so expect a
+few seconds of thinking time — which is the entire motivation for the fast path.
+
+## How it works
+
+Three tiers, each with a different accuracy and cost, described in full in
+[docs/DESIGN.md](docs/DESIGN.md).
+
+**1. Numerical simulator.** Four-state cloth dynamics (sliding, rolling,
+spinning, stationary) integrated at 1 ms, frictional ball-ball impulses with spin
+transfer and throw, cushion rebound with speed-dependent restitution, pocket
+capture, and a spatially varying cloth friction field.
+
+**2. Closed-form solver.** No integration. While a ball slides, its slip velocity
+decays along a *fixed* direction, so the friction force is constant and the path
+is exactly a parabola of known duration. Rolling is a straight line. A shot
+becomes a handful of exactly solvable segments joined at cushions, where the
+normal velocity is damped while the spin term carries through. That last detail
+is what a plain mirror-reflection approximation gets wrong: it disagreed with the
+simulator by about a metre, where this solver lands within 114 mm on direct shots
+and 225 mm across one rail.
+
+**3. Learned residual.** A small MLP predicts the vector from the closed-form
+endpoint to the simulated one. Its output head starts at zero, so training begins
+from "trust the physics" and departs only where the data insists. Its inputs
+include the closed-form solver's own conclusions — predicted endpoint, expected
+cushion count, expected pot, ghost-ball contact geometry. An ablation on the same
+architecture, epochs, seed and split puts a number on that: fed raw shot
+parameters alone it scores 179 mm on direct shots, worse than the 114 mm of the
+physics it was meant to improve; fed the solver's conclusions as well, 98 mm.
+
+## What this demonstrates
+
+Billiards is the domain; the transferable content is below.
+
+**Surrogate modelling of an expensive simulator.** The pattern — establish a
+trusted reference, find the closed-form structure inside it, learn only the
+residual, and quantify the accuracy you traded for the speed — is the same one
+used for pricing engines that are too slow for intraday risk, finite-element
+models too slow for design loops, and any Monte Carlo where the inner loop is the
+bottleneck. Almost four orders of magnitude here, and they come mostly from the
+closed form rather than from the network, which is usually how it goes.
+
+**Validating against theory rather than against yourself.** A snapshot test
+would have locked in a sign error that made the central physics wrong. Closed-form
+references, conservation laws and analytic decay rates caught it in one run. The
+same discipline applies to any numerical pipeline: solve a special case exactly
+and check against it.
+
+**Selective prediction with a free confidence signal.** The per-stratum error
+breakdown, the prediction-spread ratio and the coverage curve exist so a caller
+knows which predictions to trust — and the gating signal falls out of the physics
+baseline at no extra cost, rather than requiring a second calibrated model. A
+surrogate that is confident everywhere and accurate in half the space is more
+dangerous than the slow model it replaced.
+
+**Model selection without leaking the test set.** Epochs are chosen on a
+validation split carved from the training data; the test split is scored once,
+after training. The GBM comparison uses the identical split.
+
+**Reproducibility as a property of the code.** Dataset generation is parallel and
+seeded per sample, so the data is byte-identical on 1 core or 32. Training and
+serving construct features through one function, and a test pins the two paths
+together, because train/serve skew degrades a model quietly instead of failing
+loudly.
+
+**Porting a numerical kernel without letting it drift.** The same physics exists
+twice, in Python and in JavaScript, because the browser needed it and the
+reference had to stay the reference. A port that nobody measures is a rumour, so
+the two are pinned to each other by 35 recorded shots and compared to a
+thousandth of a millimetre in CI. The one number that made it worth doing is the
+65× speedup, which is what turns "simulate the candidate shots" from a claim
+into the thing the bot actually does inside a turn.
+
+**Deployability.** The residual model exports to ONNX and runs through ONNX
+Runtime with no PyTorch in the serving path. PyQt and OpenCV are optional extras,
+so the package installs headless. CI lints, type-checks and tests on three Python
+versions, then runs the whole data to model to benchmark to figures pipeline and
+uploads the artefacts. The playable page deploys to GitHub Pages only after the
+parity check and a run of headless games have passed.
+
+## Repository map
+
+```
+src/cueai/
+ physics/
+ ball.py four-state cloth dynamics for one ball
+ collisions.py frictional ball-ball impulses, cushions, pockets
+ simulator.py the 16-ball reference simulator
+ analytic.py the closed-form solver and its derivations
+ rack.py 8-ball rack geometry and ball identities
+ ml/
+ dataset.py parallel, per-sample-seeded shot generation
+ features.py one feature path for training and serving
+ model.py CueNet, zero-initialised residual head
+ train.py training, baseline comparison, stratified evaluation
+ infer.py predict_fast (0.60 ms) and predict (full simulation)
+ api/main.py FastAPI service
+ ui/app.py PyQt6 interactive table
+ vision/overlay.py OpenCV trajectory overlays
+
+web/ the playable table: dependency-free ES modules
+ js/
+ physics.js hand port of src/cueai/physics/, checked against it
+ rack.js the same rack geometry, ported
+ game.js eight-ball rules, fouls, group assignment
+ bot.js closed-form candidate pots, then simulated rollouts
+ aim.js ghost-ball geometry and the first-contact search
+ render.js canvas drawing: table, balls, cue, aim overlay
+ inspector.js the live cue-ball trace against 5/7·v₀
+ main.js input, the fixed-timestep loop, whose turn it is
+ test/
+ parity.mjs replays the Python reference shots, compares endpoints
+ selfplay.mjs headless bot-against-bot games, every rule branch
+ browser.mjs drives the real page in Chrome, fails on console errors
+ input.mjs plays with a real cursor and keyboard, not the test seam
+ capture.mjs records the screenshots and the clip in this README
+
+tests/
+ test_validation.py closed-form physics validation
+ test_features.py train/serve feature consistency
+ test_physics.py simulator behaviour
+ test_metrics.py the reported metrics, including the trust gate
+ test_api.py HTTP contract
+
+scripts/
+ benchmark.py writes docs/BENCHMARKS.md
+ make_figures.py writes docs/assets/*.png
+ export_parity_cases.py writes the reference shots the browser is held to
+ site_facts.py writes the measured numbers the web page quotes
+```
+
+## Limitations
+
+- Not calibrated against a real table. No measurements were taken, so the claim
+ is internal consistency with classical mechanics, not fidelity to specific
+ equipment. Cloth and cushion coefficients come from the published ranges.
+- A break under-spreads, and the size of the gap is measured rather than
+ guessed at. A rack is resolved as a chain of pairwise collisions, so
+ restitution is applied about fifteen times where the real event dissipates
+ once, and only 48% of the kinetic energy survives a 10 m/s break. Balls at the
+ back of the rack therefore leave slower than they should. This is the largest
+ known departure from reality in the simulator; the arithmetic and the fix it
+ would need are in [docs/VALIDATION.md](docs/VALIDATION.md).
+- Cue ball squirt and swerve are not modelled, so aiming advice would be
+ systematically off for heavy sidespin.
+- Multi-rail outcomes are not usefully predictable by any of the three methods
+ here, as measured above.
+- The fast path predicts *where balls stop*; it is not an aiming aid. On the layout
+ in [docs/DESIGN.md](docs/DESIGN.md) the window that pots a ball is a quarter of a
+ degree wide — `tests/test_validation.py` asserts that the ghost-ball line pots and
+ half a degree off misses — roughly ten times finer than the surrogate resolves.
+ Aiming is exact closed-form geometry anyway, and needs no model.
+- The reference simulator is straightforward Python. It is a definition of truth,
+ not a fast engine; optimisation stopped once the collision loop was no longer
+ the obvious bottleneck.
+
+## References
-> Developed a physics-informed AI simulation for billiards that combines classical
-> mechanics with regression models to predict realistic ball trajectories based on
-> spin, launch angle, velocity, cushion interactions, and table surface variations.
+The cloth and collision model follows the standard treatment in Ron Shepard's
+*Amateur Physics for the Amateur Pool Player*, Wayland Marlow's *The Physics of
+Pocket Billiards*, and David Alciatore's technical proofs; the four-state
+formulation matches the approach taken by
+[pooltool](https://github.com/ekiefl/pooltool). Coefficients are from the
+published ranges in those sources.
## License
-MIT
+MIT — see [LICENSE](LICENSE).
diff --git a/docs/BENCHMARKS.md b/docs/BENCHMARKS.md
new file mode 100644
index 0000000..770a460
--- /dev/null
+++ b/docs/BENCHMARKS.md
@@ -0,0 +1,69 @@
+# Benchmarks
+
+Regenerate with `python scripts/benchmark.py`. Single CPU core, no GPU.
+
+- Python 3.12.3 on Linux-6.12.94+-x86_64-with-glibc2.39
+- Processor: x86_64
+
+## Latency per shot
+
+| Method | Mean | Median | p95 | Speedup vs full rack |
+|---|---:|---:|---:|---:|
+| Numerical simulator, 16-ball rack | 4.6 s | 4.6 s | 4.6 s | — |
+| Numerical simulator, cue + object ball | 337.9 ms | 337.8 ms | 338.7 ms | 14x |
+| Closed-form solver (no ML) | 0.266 ms | 0.264 ms | 0.275 ms | 17,290x |
+| Gradient boosting on the same features | 2.6 ms | 2.6 ms | 2.8 ms | 1,744x |
+| Closed form + CueNet residual (ONNX Runtime) | 0.602 ms | 0.572 ms | 0.590 ms | 7,638x |
+| CueNet forward pass only, batch of 1024 (per shot) | 1.6 us | 0.8 us | 4.7 us | — |
+
+## Accuracy on held-out shots
+
+4,000 test shots from a 20,000 shot dataset. Error is the distance between predicted and simulated resting position, averaged over the cue ball and the object ball.
+
+| Model | Mean error | p95 | Cue ball | Object ball | R² |
+|---|---:|---:|---:|---:|---:|
+| Closed form (no fitting) | 494 mm | 1778 mm | 730 mm | 258 mm | 0.013 |
+| Gradient boosting on raw features | 382 mm | 1163 mm | 516 mm | 249 mm | 0.521 |
+| Closed form + CueNet residual | 376 mm | 1293 mm | 516 mm | 236 mm | 0.431 |
+
+## Where prediction stops working
+
+Resting position is a smooth function of the shot until the ball starts ricocheting between cushions. Mean error in mm by cushion contacts:
+
+| Cushion contacts | Shots | Closed form | Gradient boosting | CueNet residual |
+|---|---:|---:|---:|---:|
+| 0 | 414 | 114 | 162 | 98 |
+| 1 | 650 | 225 | 222 | 187 |
+| 2 | 876 | 371 | 289 | 281 |
+| 3+ | 2060 | 708 | 517 | 532 |
+
+## Reading the object-ball number honestly
+
+Only about a third of sampled shots actually hit the object ball. For the rest it stays where it started, which the closed-form baseline predicts exactly, so those rows contribute almost no object-ball error to any model. Split out, in mm:
+
+| Ball-ball contact | Shots | Share | Closed form cue / object | Gradient boosting cue / object | CueNet cue / object |
+|---|---:|---:|---:|---:|---:|
+| no | 2705 | 67.6% | 585 / 0 | 457 / 77 | 411 / 16 |
+| yes | 1295 | 32.4% | 1033 / 798 | 638 / 610 | 736 / 696 |
+
+## Feature ablation
+
+Same architecture, epochs, seed and split; the ablated model sees only raw shot parameters, with none of the closed-form solver's conclusions. This is the measurement behind the claim that the features carry this result:
+
+| CueNet inputs | All shots | Direct shots (no cushion) |
+|---|---:|---:|
+| Raw shot parameters only | 469 mm | 179 mm |
+| Plus closed-form output | 376 mm | 98 mm |
+| _closed form alone, for reference_ | 494 mm | 114 mm |
+
+## Knowing which predictions to trust
+
+The cushion-contact breakdown above slices by what the simulator did, which is only knowable after paying for the simulation. This table slices by the cushion count the closed-form solver *expects*, which is already computed as part of making the prediction and therefore costs nothing extra. Answering only the shots below a threshold trades coverage for accuracy:
+
+| Expected cushions | Shots answered | Coverage | Closed form | Gradient boosting | CueNet residual |
+|---|---:|---:|---:|---:|---:|
+| ≤ 0 | 390 | 9.8% | 113 mm | 160 mm | 100 mm |
+| ≤ 1 | 1098 | 27.5% | 240 mm | 221 mm | 189 mm |
+| ≤ 2 | 2002 | 50.0% | 334 mm | 265 mm | 251 mm |
+| ≤ 3 | 2748 | 68.7% | 401 mm | 314 mm | 303 mm |
+| no gate | 4000 | 100.0% | 494 mm | 382 mm | 376 mm |
diff --git a/docs/DESIGN.md b/docs/DESIGN.md
new file mode 100644
index 0000000..1a0eef1
--- /dev/null
+++ b/docs/DESIGN.md
@@ -0,0 +1,250 @@
+# Design
+
+## The problem this solves
+
+Simulating a billiards shot accurately is slow. Sixteen bodies, each in one of
+four motion states, with collisions resolved every millisecond, costs about
+**4.6 seconds per rack shot** in this implementation. Anything that needs many
+shot outcomes — searching for a good shot, estimating the odds of a break,
+training a policy — cannot afford that.
+
+The usual answer is to replace the simulator with a learned model. Done naively
+that throws away physics that is already exact, and it produces a model whose
+errors have no structure. The approach here keeps the physics and learns only
+what the physics leaves out.
+
+```
+shot parameters (speed, angle, tip offset, cue position, cloth μ, cushion e)
+ |
+ +-----------------------------+
+ | |
+ v v
+ closed-form solver CueNet residual
+ (exact segments, 0.27 ms) (what the closed form misses)
+ | |
+ +--------------+--------------+
+ v
+ predicted resting positions 0.60 ms total
+ :
+ : compared against
+ v
+ numerical simulator 4.6 s, ground truth
+```
+
+## Three tiers, on purpose
+
+**Tier 1: the numerical simulator** (`cueai.physics.simulator`). Explicit Euler
+integration at 1 ms, four-state cloth dynamics, frictional ball-ball impulses
+with spin transfer, cushion rebound, pocket capture. This is the definition of
+truth for everything else, and it is validated against closed-form results in
+[VALIDATION.md](VALIDATION.md).
+
+**Tier 2: the closed-form solver** (`cueai.physics.analytic`). No integration at
+all. The observation that makes this possible: while a ball slides, the slip
+velocity `u` decays along a *fixed direction*, so the friction force is constant
+and the path over that phase is exactly a parabola of known duration
+`|u|/(3.5 μ_s g)`. Rolling is a straight line of length `v²/(2 μ_r g)`. A shot is
+therefore a short sequence of exactly solvable segments joined at cushion
+contacts, where the normal velocity is damped and the spin term `u - v` carries
+through — which is why a ball leaves a rail sliding rather than rolling. Getting
+that detail right matters: a plain mirror-reflection approximation, which assumes
+the ball leaves the rail rolling, disagreed with the simulator by about a metre,
+while this solver lands within 114 mm on direct shots and 225 mm across one rail.
+
+**Tier 3: the learned residual** (`cueai.ml`). A small MLP predicts the vector
+from the closed-form endpoint to the simulated endpoint. Its head is initialised
+to zero, so training starts from "trust the physics exactly" and moves away only
+where the data insists.
+
+## Why the features matter more than the architecture
+
+The first working version of this model made predictions *worse* than the
+closed-form solver on the easiest shots. The network was being fed raw shot
+parameters, so to know where the closed-form model went wrong it first had to
+rediscover cushion reflection geometry from a speed and an angle. It spent its
+capacity on the chaotic shots that dominate the loss and degraded the clean ones.
+
+The fix was to hand the network what the physics already knew: the closed-form
+endpoint, the cushion count it expects, whether it expects a pot, and the
+ghost-ball geometry deciding whether the object ball is contacted at all.
+
+Every training run re-measures this as an ablation, on the same architecture,
+epochs, seed and split, so the claim is not a story about an earlier commit:
+
+| CueNet inputs | All shots | Direct shots |
+|---|---:|---:|
+| Raw shot parameters only | 469 mm | 179 mm |
+| Plus the closed-form solver's output | 376 mm | 98 mm |
+| _closed form alone, for reference_ | 494 mm | 114 mm |
+
+Without the physics features the network is *worse than the physics* on the shots
+the physics nearly solves — 179 mm against 114 mm — while still beating it on
+average, which is exactly the failure that a single headline number conceals. With
+them it is better on both. The predicted cushion count does double duty: it also
+lets the model recognise a chaotic shot and hedge instead of guessing, and it is
+the gate in the next section.
+
+## Where it stops working
+
+A shot's resting position is a smooth function of its inputs right up to the
+point where the ball starts ricocheting. After a few cushion contacts, a
+millimetre of difference in the cue position moves the final position by a table
+length, and no model recovers that.
+
+
+
+So the results are reported by cushion contact count rather than as a single
+average. The residual model is the best of the three for shots with 0, 1 or 2
+cushion contacts; past that, plain gradient boosting edges ahead by hedging
+harder toward the middle of the table. The reported spread ratio makes that
+visible: CueNet reproduces 97% of the true spread in the predictable buckets and
+86% in the chaotic one.
+
+## Reporting confidence without a second model
+
+Slicing by the simulator's cushion count is an autopsy: you only have that number
+once you have paid for the simulation. The closed-form solver, though, reports the
+cushion count it *expects* while computing its answer, so that signal is free —
+and it is a good enough proxy for "is this outcome chaotic" to gate on:
+
+| Answer only when the solver expects | Coverage | CueNet error |
+|---|---:|---:|
+| No cushion | 9.8% | 100 mm |
+| At most one | 27.5% | 189 mm |
+| At most two | 50.0% | 251 mm |
+| Anything | 100% | 376 mm |
+
+
+
+This is what makes the fast path deployable rather than merely fast: a caller can
+choose an error budget, take the coverage that comes with it, and route the
+remainder to the simulator. The alternative — a single headline error over a
+distribution that mixes 100 mm shots with 700 mm shots — gives a caller no way to
+act. The residual model happens to be the most accurate of the three at every
+coverage level, which is a more robust claim than its 1.6% edge on the mean.
+
+## What the surrogate is and is not good enough for
+
+Worth being precise, because "0.60 ms billiards predictor" invites the wrong
+conclusion.
+
+**Not accurate enough to aim with.** On a representative layout — cue ball at
+(0.60, 0.35), object ball at (1.70, 0.75), far corner pocket — the aim window
+that actually pots the ball is **0.15° to 0.25° wide**, at stroke speeds of 1.5,
+2.5 and 4 m/s. Ranking thousands of candidate aim lines by the surrogate's
+predicted object endpoint does concentrate them near the correct angle, but the
+top-ranked candidates spread over roughly two degrees, so the ranking never
+resolves a window an order of magnitude narrower than that.
+
+This is less of a loss than it looks, because aiming needs no learned model at
+all. The contact geometry is exact and closed-form: the ghost-ball point is the
+object ball's centre displaced one diameter back along the line to the pocket, and
+shooting at it pots the ball at every speed tested — which
+`tests/test_validation.py` asserts, along with half a degree off it missing. The
+window is not symmetric
+about it, though — it extends about 0.1° to the thin side and barely at all to the
+thick side, which is collision-induced throw, the object ball being dragged off
+the geometric line by ball-ball friction. That the simulator reproduces an effect
+real players compensate for by feel is a better argument for the cloth model than
+any of the aggregate error numbers.
+
+**Weakest exactly where the baseline is weakest.** Two thirds of sampled shots
+never reach the object ball, and for those the baseline's "it stays put" is exactly
+right, so the reported object-ball error of 236 mm is mostly an average over free
+zeros. On the third of shots that do involve a collision, the numbers are 736 mm
+for the cue ball and 696 mm for the object ball — and plain gradient boosting beats
+the residual model there, 638 mm and 610 mm.
+
+That is a structural consequence of the design rather than a tuning problem. The
+closed-form solver has no ball-ball contact model, so on a collision it can be a
+metre wrong, and the residual is then asked to undo a large error rather than
+refine a small one. A model predicting endpoints directly has no such anchor to
+fight. The corollary is that the highest-value next change is not a bigger network
+or more data: it is giving the closed-form solver a ghost-ball collision so the
+baseline it hands over is worth correcting.
+
+The residual formulation wins the other side of that trade. On shots where nothing
+happens it introduces 16 mm of spurious object-ball motion against gradient
+boosting's 77 mm, because "predict no correction" is where it starts rather than
+something it has to infer.
+
+So the honest summary: a good surrogate for the *distribution* of low-cushion
+outcomes and for screening large candidate sets, an exact tool for contact
+geometry, weaker than a plain regressor once a collision is involved, and not a
+substitute for the simulator when a specific multi-rail outcome matters.
+
+## The fourth tier: the same physics, in a browser
+
+Nothing above can be watched, and a simulator that cannot be watched is taken on
+trust. `web/` is a playable eight-ball table running the same model, which turns
+every claim in this document into something a reader can check by shooting.
+
+**Why a hand port rather than Pyodide or WebAssembly.** Shipping the Python
+would have kept one implementation, which is the obvious argument for it. It
+also ships a multi-megabyte runtime for a physics loop of a few hundred lines,
+and it makes the page's responsiveness a property of someone else's interpreter.
+The loop is small, it is the part of the system whose behaviour is best pinned
+down, and it needs to run tens of thousands of steps inside a turn. So it was
+ported by hand, and the cost of that decision — two implementations that can
+drift — is paid for directly rather than hoped away.
+
+**How the port is held to the reference.** `scripts/export_parity_cases.py`
+records 35 shots from the Python simulator: draw, follow, english off two rails,
+thin cuts, clusters, and full sixteen-ball breaks, chosen so that everything the
+model does appears in at least one of them. `web/test/parity.mjs` replays each
+one in Node and compares every ball's resting position. The two agree to
+1.2 × 10⁻³ mm across 156 s of simulated table time, which is not a tolerance so
+much as a demonstration that the same arithmetic is being done in the same
+order. It runs in CI, and the page does not deploy if it fails.
+
+**The opponent is a search, which is the point.** Aiming needs no model: the
+ghost-ball construction is exact, and a test asserts it pots while half a degree
+either side misses. What is actually hard is *choosing* which exact shot to
+play, and that is a question about futures rather than geometry. So `web/js/bot.js`
+enumerates every ball-and-pocket pair in closed form, discards what is blocked
+or cut thinner than 78°, orders the survivors by a cheap pottability prior, and
+then spends its entire budget simulating them — scoring each rollout by what it
+leaves behind, not only by whether it drops. The strength setting is a rollout
+budget, and the panel reports it in rollouts and table time rather than as an
+adjective, because that is the honest unit.
+
+**The search does not use the surrogate, and that is a result.** A surrogate
+earns its error when you must screen far more candidates than you can afford to
+simulate. Inside one turn there are a few dozen candidates and the ported
+simulator runs 65× faster than the Python reference, so the exact answer is
+affordable and the approximate one would only add error. The place the surrogate
+would earn its keep is the case this game does not present: sweeping a
+continuous space of speed, angle and spin per candidate, where the count is
+bounded by the budget rather than by the geometry.
+
+**Two invariants that only a long run can check.** Twenty headless
+bot-against-bot games (`web/test/selfplay.mjs`) walk every branch of the rules,
+and on each of roughly 800,000 physics steps assert that no two balls share
+space and nothing is inside a cushion. The worst overlap seen is 0.47 mm on a
+57 mm ball. This is the check that distinguishes a solver that converges from
+one that merely looks like it does, and no single shot would ever reveal it.
+
+## Choices a reviewer might question
+
+**Explicit Euler at 1 ms rather than something higher order.** The dynamics are
+piecewise smooth with impulsive events at contacts, so the local error is
+dominated by event handling, not by integrator order. The validation suite bounds
+the resulting drift at under 1% against closed form. A symplectic or adaptive
+integrator would be the right call if the cloth model were stiffer.
+
+**Positional correction on overlap.** Overlapping balls are pushed apart before
+the impulse is applied, which is a projection rather than a physical force. It is
+what keeps a packed rack stable; the alternative is sub-stepping to the exact
+contact time, which costs more than the accuracy is worth here.
+
+**Endpoints as the learning target.** Resting position is what shot selection
+needs. Predicting whole trajectories would be a sequence-model problem, and the
+chaos analysis above suggests the horizon over which that is worth attempting is
+short.
+
+**Dataset size.** 20,000 shots, which is about fifteen minutes of generation on
+eight cores. No learning curve was measured, so this is a practical choice rather
+than a demonstrated sufficiency — the honest expectation is that more data would
+help the collision cases, where the model is fighting the baseline, and do very
+little for the multi-rail cases, where the target is chaotic. Generation is
+parallel and seeded per sample, so the dataset is identical on 1 core or 32.
diff --git a/docs/VALIDATION.md b/docs/VALIDATION.md
new file mode 100644
index 0000000..21af412
--- /dev/null
+++ b/docs/VALIDATION.md
@@ -0,0 +1,196 @@
+# Validation
+
+The simulator is checked against closed-form solutions and conservation laws
+rather than against its own previous output, so a regression means the physics is
+wrong and not merely different. Everything below is asserted in
+`tests/test_validation.py` and runs in about ten seconds.
+
+```bash
+pytest tests/test_validation.py -v
+```
+
+
+
+## The cloth model
+
+A ball on cloth is in one of four states, decided by the velocity of its contact
+point, `u = v + ω × (-Rẑ)`:
+
+| State | Condition | Dynamics |
+|---|---|---|
+| Sliding | `u ≠ 0` | friction `μ_s g` opposite `u`, with the matching torque |
+| Rolling | `u = 0` | rolling resistance `μ_r g` opposite `v` |
+| Spinning | `v = 0`, `ω_z ≠ 0` | `dω_z/dt = -5 μ_sp g / 2R` |
+| Stationary | both zero | at rest |
+
+While sliding, the slip velocity decays 3.5x faster than the centre of mass,
+because friction both slows the ball and spins it up:
+
+```
+|du/dt| = μ_s g (1 + mR²/I) = μ_s g (1 + 5/2)
+```
+
+## What is checked
+
+| Property | Reference | Tolerance | Measured |
+|---|---|---|---|
+| Speed when rolling begins | `5/7 v₀` | 1% | 0.3% |
+| Sliding distance | `12v₀² / (49 μ_s g)` | 2% | 0.8% |
+| Total stopping distance | slide + `v²/(2 μ_r g)` | 1% | 0.6% |
+| Natural roll tip offset | no sliding phase at `f = 0.4` | exact | `u = 0` |
+| Vertical spin decay | `5 μ_sp g / 2R` | 1e-6 | exact |
+| Cushion rebound speed | `e · v_approach` | 1e-9 | exact |
+| Momentum in a frictional collision | conserved | 1e-12 | exact |
+| Energy in any contact | never increases | — | holds |
+| Sidespin on a rolling ball | no lateral deflection | 1e-9 m | exact |
+| Draw / stun / follow | ordered, draw ends behind contact | — | holds |
+| Ghost-ball aim line | pots the ball; ±0.5° misses | — | holds |
+| Timestep convergence, full-ball contact | 2 ms vs 1 ms | 10 mm | 3 mm |
+| Timestep convergence, thin cut | 2 ms vs 1 ms | 60 mm | 46 mm |
+| Contacts a racked triangle presents | 30 | exact | 30 |
+| Balls set moving by a full-power break | all 15 | — | holds |
+| Break spread against cue speed | strictly increasing | — | holds |
+
+## The bug none of the above could see
+
+Every property in that table concerns one ball, or two. A defect that lives in
+the relationship between fifteen of them passes all of it, and one did.
+
+Two balls were treated as touching when the gap between their surfaces was at
+most `1e-4 m`. The rack was built with a `1e-4 m` clearance. So all thirty
+contacts in the triangle sat exactly on the threshold that decides whether a
+contact exists, and which side each one landed on was decided by whether
+`hypot` rounded up or down for that particular pair. Sixteen registered.
+Fourteen did not.
+
+The consequences were visible only in aggregate. A break propagated through a
+contact graph with holes in it, so balls in the middle of the rack came out of
+a full-power break having barely moved, and the table opened up *less* the
+harder it was struck. Both configurations, 20 racks at each speed:
+
+| Cue speed | 3 m/s | 5 m/s | 6.5 m/s | 8.2 m/s | 10 m/s |
+|---|---:|---:|---:|---:|---:|
+| Mean pair separation, before | 0.348 m | 0.610 m | 0.575 m | **0.478 m** | 0.577 m |
+| Mean pair separation, after | 0.320 m | 0.568 m | 0.577 m | 0.637 m | 0.731 m |
+| Mean distance travelled, before | 0.162 m | 0.321 m | 0.309 m | 0.310 m | 0.387 m |
+| Mean distance travelled, after | 0.135 m | 0.303 m | 0.370 m | 0.474 m | 0.549 m |
+
+Before the fix, hitting the rack half again as hard as 5 m/s left it *tighter*,
+and past 5 m/s the balls stopped travelling any further at all. After it, both
+rise with every increase in cue speed.
+
+It was found by measuring that relationship and getting the sign wrong, not by
+a test going red — and not by the browser parity harness either, which
+reproduced the broken contact graph to eleven decimal places because it was a
+faithful port of it. Two implementations agreeing is evidence about the port,
+and about nothing else.
+
+The tolerance now has a name, `CONTACT_BAND`, the broad and narrow phases use
+the same one, and the balls are racked touching so that nothing sits on the
+boundary. Tests assert all three, including that the band stays far from both
+the floating-point noise that would swallow it and the physical scale that
+would make it meaningless — and, now, that a harder break opens the table
+further, which is the property nobody thought to assert because it is too
+obvious to state.
+
+## The largest known departure from reality
+
+Separately from that bug, and not fixed by it: resolving a rack's contacts
+pairwise, in sequence, dissipates far more energy than a real break does.
+
+| Cue speed | 3 m/s | 5 m/s | 6.5 m/s | 8.2 m/s | 10 m/s |
+|---|---:|---:|---:|---:|---:|
+| Kinetic energy surviving the impact | 74% | 65% | 59% | 53% | 48% |
+
+Measured either side of the largest single-step energy drop in the shot, which
+is the moment the rack is struck; before and after the contact-band fix these
+numbers agree to within a point, so this is not a consequence of it.
+
+The mechanism is arithmetic rather than mysterious. An equal-mass head-on
+collision with restitution `e` keeps `(1 + e²)/2` of the kinetic energy, about
+95% at `e = 0.95`. A rack is fifteen balls already touching, so the impulse is
+resolved as a chain of roughly fifteen such collisions and the survival factor
+compounds to `0.95¹⁵ ≈ 0.46`. The real event is a single stress wave through
+bodies in continuous contact, which dissipates once rather than fifteen times.
+
+What it costs: the speed the impulse hands to each ball falls off geometrically
+down the chain, so the balls at the back of the rack leave with a few
+centimetres per second and finish near where they started. A break here opens
+the table, but it opens it less than a real one, and the difference grows with
+cue speed. Fixing it properly means treating a simultaneous contact set as one
+event — applying material restitution at the leading edge of the propagation
+rather than at every pair along it — which is a change to the solver rather
+than to a coefficient, and is not attempted here.
+
+## Is the reference converged?
+
+Worth asking, because everything downstream is measured against it. Training
+labels are generated at a 2 ms step, so if that step were too coarse the reported
+model errors would partly be measuring the integrator.
+
+Halving the step moves the resting position by about 3 mm on a full-ball contact
+and 46 mm on a thin cut, the thin cut being the worst case because a small change
+in contact geometry is amplified into a large change in direction. Quartering it,
+to 0.5 ms, moves free-ball endpoints by a median of 0.8 mm and at most 19 mm.
+Against the 100 mm error the surrogate achieves on its best shots, discretisation
+is a real but sub-dominant term — and it is bounded by a test rather than assumed.
+
+The tip offset result is worth spelling out because it is easy to get wrong. A
+horizontal impulse applied a distance `d = f·R` off centre gives `Δv = J/m` and
+`Δω = J d / I`, so with `I = (2/5)mR²`:
+
+```
+ω = 2.5 · f · v / R
+```
+
+Pure rolling requires `ω = v/R`, so `f = 0.4`: a tip `2R/5` above centre, `7R/5`
+above the cloth, launches the ball already rolling. Offsets beyond `|f| = 0.5`
+are past the practical miscue limit, which is where the dataset sampling stops.
+
+
+
+Draw deserves a note as well. A ball struck with backspin on an open table does
+**not** come back: friction removes the backspin before it removes the forward
+velocity, and it ends up rolling forward. Draw only works because a full-ball
+collision takes the forward velocity away while the backspin survives. The test
+therefore asserts the ordering across a contact, not the behaviour of a free
+ball, and the figure above shows the cue ball finishing 0.18 m up-table, well
+behind the 1.14 m contact point.
+
+## What is not modelled
+
+Stated plainly, because a validation document that only lists successes is not
+worth much:
+
+- **Cue ball squirt and swerve.** A real off-centre hit deflects the cue ball a
+ degree or two off the aim line and then curves it. Neither is modelled, so
+ aiming corrections from this simulator would be systematically wrong for
+ heavy english.
+- **Massé.** Cue elevation adds a little vertical-axis spin and nothing else.
+- **Cushion geometry.** Rails are treated as flat vertical planes at the ball
+ centre height, with a speed-dependent restitution as a stand-in for cushion
+ compliance. Real cushions contact above centre and lift the ball slightly.
+- **Pocket geometry.** A pocket is a capture radius, not a jaw. Balls that would
+ rattle out in reality drop here.
+- **Cloth inhomogeneity** is a smooth analytic noise field, not a measured one.
+- **Ball-ball throw** uses a velocity-dependent friction coefficient of the
+ usual form `μ_b ≈ a + b·exp(-c|v_rel|)`, with coefficients that are plausible
+ rather than fitted to measurements.
+- **Ball-ball restitution is treated as speed-dependent**, falling from 0.95 by
+ 0.02 for every m/s of approach speed above 2. The cushion's speed dependence
+ is measured in the literature; this one is not, and it is the reason a break
+ keeps 48% of its energy at 10 m/s rather than 75%. Removing it was tried:
+ energy retention becomes flat in speed, as it should be if restitution is,
+ and mean pair separation after a 10 m/s break rises from 0.73 m to 0.79 m —
+ but the *median* ball moves less, not more, because a chain of near-elastic
+ collisions transmits more completely to the ball at the end of it and leaves
+ the middle of the rack even quieter. It is left in place because changing it
+ would be tuning an uncalibrated parameter to no measured benefit, but it is a
+ choice rather than a result.
+- **Simultaneous multi-ball contact**, as described above: a rack is resolved as
+ a sequence of pairwise collisions, which over-counts the dissipation and
+ under-spreads a break.
+
+None of these are calibrated against a real table, because no measurements were
+taken. The claim this project makes is internal consistency with classical
+mechanics, not fidelity to a specific piece of equipment.
diff --git a/docs/assets/accuracy.png b/docs/assets/accuracy.png
new file mode 100644
index 0000000..965e3d7
Binary files /dev/null and b/docs/assets/accuracy.png differ
diff --git a/docs/assets/break_shot.png b/docs/assets/break_shot.png
new file mode 100644
index 0000000..ab5d33e
Binary files /dev/null and b/docs/assets/break_shot.png differ
diff --git a/docs/assets/latency.png b/docs/assets/latency.png
new file mode 100644
index 0000000..18be704
Binary files /dev/null and b/docs/assets/latency.png differ
diff --git a/docs/assets/reliability.png b/docs/assets/reliability.png
new file mode 100644
index 0000000..d9f9216
Binary files /dev/null and b/docs/assets/reliability.png differ
diff --git a/docs/assets/spin_control.png b/docs/assets/spin_control.png
new file mode 100644
index 0000000..6b2702e
Binary files /dev/null and b/docs/assets/spin_control.png differ
diff --git a/docs/assets/validation.png b/docs/assets/validation.png
new file mode 100644
index 0000000..cfa0942
Binary files /dev/null and b/docs/assets/validation.png differ
diff --git a/docs/assets/web_bot.png b/docs/assets/web_bot.png
new file mode 100644
index 0000000..cb4a04a
Binary files /dev/null and b/docs/assets/web_bot.png differ
diff --git a/docs/assets/web_demo.gif b/docs/assets/web_demo.gif
new file mode 100644
index 0000000..1502fe0
Binary files /dev/null and b/docs/assets/web_demo.gif differ
diff --git a/docs/assets/web_demo.mp4 b/docs/assets/web_demo.mp4
new file mode 100644
index 0000000..68a01f1
Binary files /dev/null and b/docs/assets/web_demo.mp4 differ
diff --git a/docs/assets/web_game.png b/docs/assets/web_game.png
new file mode 100644
index 0000000..1a89014
Binary files /dev/null and b/docs/assets/web_game.png differ
diff --git a/docs/assets/web_inspector.png b/docs/assets/web_inspector.png
new file mode 100644
index 0000000..3e94777
Binary files /dev/null and b/docs/assets/web_inspector.png differ
diff --git a/models/latency.json b/models/latency.json
new file mode 100644
index 0000000..3302321
--- /dev/null
+++ b/models/latency.json
@@ -0,0 +1,43 @@
+{
+ "simulator_full_rack": {
+ "mean_ms": 4595.944689000437,
+ "median_ms": 4595.326960999955,
+ "p95_ms": 4598.883345201102,
+ "n": 3
+ },
+ "simulator_two_ball": {
+ "mean_ms": 337.9245412999808,
+ "median_ms": 337.8439725001954,
+ "p95_ms": 338.7264999002582,
+ "n": 10,
+ "speedup_vs_full_rack": 13.600505815055756
+ },
+ "closed_form": {
+ "mean_ms": 0.2658112820063252,
+ "median_ms": 0.26385650016891304,
+ "p95_ms": 0.2752571504061052,
+ "n": 500,
+ "speedup_vs_full_rack": 17290.254402712195
+ },
+ "gradient_boosting": {
+ "mean_ms": 2.634616043946153,
+ "median_ms": 2.601596499516745,
+ "p95_ms": 2.776304000781238,
+ "n": 500,
+ "speedup_vs_full_rack": 1744.4457227689952
+ },
+ "surrogate_onnx": {
+ "mean_ms": 0.6016930600089836,
+ "median_ms": 0.5721709994759294,
+ "p95_ms": 0.5901895490751485,
+ "n": 500,
+ "speedup_vs_full_rack": 7638.354161724614
+ },
+ "cuenet_batch1024": {
+ "mean_ms": 0.0015760251171315076,
+ "median_ms": 0.0007746210926029562,
+ "p95_ms": 0.004728075389692776,
+ "n": 25,
+ "speedup_vs_full_rack": 2916162.0833590687
+ }
+}
diff --git a/models/metrics.json b/models/metrics.json
index 182a377..9d66bce 100644
--- a/models/metrics.json
+++ b/models/metrics.json
@@ -1,10 +1,148 @@
{
- "torch_mae_m": 0.007152832578867674,
- "physics_only_mae_m": 0.006510819774121046,
- "improvement_pct": -9.86070613255915,
- "epochs": 12,
- "n_samples": 600,
- "history_last": 2.6276023224151383e-05,
- "onnx": "models/cuenet.onnx",
- "checkpoint": "models/cuenet.pt"
-}
\ No newline at end of file
+ "n_samples": 20000,
+ "n_test": 4000,
+ "epochs": 300,
+ "hidden": 256,
+ "final_train_loss": 0.05557978451251984,
+ "models": {
+ "analytic": {
+ "mae_mm": 308.66009521484375,
+ "euclidean_mm": 494.02587890625,
+ "p95_mm": 1777.5478515625,
+ "cue_mm": 729.7384033203125,
+ "obj_mm": 258.3133544921875,
+ "r2": 0.013224394297231212
+ },
+ "gbm": {
+ "mae_mm": 243.47491761920614,
+ "euclidean_mm": 382.44517827535526,
+ "p95_mm": 1163.1076422861256,
+ "cue_mm": 515.7021449395502,
+ "obj_mm": 249.18821161116028,
+ "r2": 0.5207620486955216
+ },
+ "cuenet": {
+ "mae_mm": 237.5008544921875,
+ "euclidean_mm": 376.185791015625,
+ "p95_mm": 1292.6614990234375,
+ "cue_mm": 516.0587158203125,
+ "obj_mm": 236.31285095214844,
+ "r2": 0.4314108586578608
+ }
+ },
+ "by_cushion_contacts": [
+ {
+ "cushion_contacts": "0",
+ "n": 414,
+ "analytic": 113.5447998046875,
+ "gbm": 161.661007928057,
+ "cuenet": 98.3550033569336,
+ "cuenet_spread_ratio": 0.9720291495323181
+ },
+ {
+ "cushion_contacts": "1",
+ "n": 650,
+ "analytic": 224.95225524902344,
+ "gbm": 222.43481668532624,
+ "cuenet": 186.57786560058594,
+ "cuenet_spread_ratio": 0.9690799713134766
+ },
+ {
+ "cushion_contacts": "2",
+ "n": 876,
+ "analytic": 370.912109375,
+ "gbm": 289.2407014788066,
+ "cuenet": 281.2088928222656,
+ "cuenet_spread_ratio": 0.952953577041626
+ },
+ {
+ "cushion_contacts": "3+",
+ "n": 2060,
+ "analytic": 707.74658203125,
+ "gbm": 516.9395973195673,
+ "cuenet": 532.2376708984375,
+ "cuenet_spread_ratio": 0.8552145957946777
+ }
+ ],
+ "by_object_ball_contact": [
+ {
+ "object_ball_contact": "no",
+ "n": 2705,
+ "share_pct": 67.6,
+ "analytic_cue_mm": 584.77587890625,
+ "analytic_obj_mm": 1.0870197002077475e-05,
+ "gbm_cue_mm": 457.21557076675373,
+ "gbm_obj_mm": 76.60472891374283,
+ "cuenet_cue_mm": 410.94439697265625,
+ "cuenet_obj_mm": 16.42076873779297
+ },
+ {
+ "object_ball_contact": "yes",
+ "n": 1295,
+ "share_pct": 32.4,
+ "analytic_cue_mm": 1032.5364990234375,
+ "analytic_obj_mm": 797.8789672851562,
+ "gbm_cue_mm": 637.8690817251982,
+ "gbm_obj_mm": 609.6811233459201,
+ "cuenet_cue_mm": 735.621826171875,
+ "cuenet_obj_mm": 695.6240844726562
+ }
+ ],
+ "risk_coverage": [
+ {
+ "expected_cushions_at_most": "0",
+ "n": 390,
+ "coverage_pct": 9.8,
+ "analytic": 113.08963775634766,
+ "gbm": 159.65888517893276,
+ "cuenet": 99.91382598876953
+ },
+ {
+ "expected_cushions_at_most": "1",
+ "n": 1098,
+ "coverage_pct": 27.5,
+ "analytic": 240.04368591308594,
+ "gbm": 220.96904733242332,
+ "cuenet": 189.33424377441406
+ },
+ {
+ "expected_cushions_at_most": "2",
+ "n": 2002,
+ "coverage_pct": 50.0,
+ "analytic": 333.79351806640625,
+ "gbm": 265.3830192910306,
+ "cuenet": 250.65756225585938
+ },
+ {
+ "expected_cushions_at_most": "3",
+ "n": 2748,
+ "coverage_pct": 68.7,
+ "analytic": 401.277587890625,
+ "gbm": 314.3183330140579,
+ "cuenet": 303.30987548828125
+ },
+ {
+ "expected_cushions_at_most": "all",
+ "n": 4000,
+ "coverage_pct": 100.0,
+ "analytic": 494.02587890625,
+ "gbm": 382.44517827535526,
+ "cuenet": 376.185791015625
+ }
+ ],
+ "feature_ablation": {
+ "note": "identical architecture, epochs, seed and split; the ablated model sees only raw shot parameters, not the closed-form solver's output",
+ "shot_features_only_mm": 469.4621276855469,
+ "with_closed_form_features_mm": 376.185791015625,
+ "shot_features_only_direct_mm": 179.4994659423828,
+ "with_closed_form_features_direct_mm": 98.3550033569336,
+ "analytic_direct_mm": 113.5447998046875
+ },
+ "error_reduction_vs_analytic_pct": 23.9,
+ "error_reduction_vs_gbm_pct": 1.6,
+ "environment": {
+ "python": "3.12.3",
+ "torch": "2.13.0+cpu",
+ "platform": "Linux-6.12.94+-x86_64-with-glibc2.39"
+ }
+}
diff --git a/pyproject.toml b/pyproject.toml
index b5ffa9b..3a2ad25 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,28 +4,37 @@ build-backend = "setuptools.build_meta"
[project]
name = "cueai"
-version = "0.1.0"
-description = "Physics-informed AI billiards simulation platform"
+version = "0.3.0"
+description = "Physics-informed billiards simulation with a learned fast surrogate"
readme = "README.md"
-requires-python = ">=3.9"
+requires-python = ">=3.10"
+license = { text = "MIT" }
dependencies = [
"numpy>=1.26",
"pandas>=2.1",
- "scipy>=1.11",
"scikit-learn>=1.4",
+ "joblib>=1.3",
"torch>=2.1",
"onnx>=1.15",
"onnxruntime>=1.16",
- "opencv-python-headless>=4.8",
- "PyQt6>=6.6",
"fastapi>=0.109",
"uvicorn[standard]>=0.27",
"pydantic>=2.5",
- "httpx>=0.26",
- "matplotlib>=3.8",
"tqdm>=4.66",
]
+[project.optional-dependencies]
+# Desktop UI and OpenCV overlays are not needed to train, serve or test.
+ui = ["PyQt6>=6.6"]
+vision = ["opencv-python-headless>=4.8"]
+dev = [
+ "pytest>=7.4",
+ "ruff>=0.5",
+ "mypy>=1.8",
+ "matplotlib>=3.8",
+ "httpx>=0.26",
+]
+
[project.scripts]
cueai-ui = "cueai.ui.app:main"
cueai-train = "cueai.ml.train:main"
@@ -37,3 +46,31 @@ where = ["src"]
[tool.pytest.ini_options]
testpaths = ["tests"]
pythonpath = ["src"]
+
+[tool.ruff]
+line-length = 100
+target-version = "py310"
+src = ["src", "tests", "scripts"]
+
+[tool.ruff.lint]
+select = ["E", "F", "W", "I", "UP", "B", "C4", "SIM", "RUF"]
+ignore = [
+ "E741", # physics notation: I for inertia, e for restitution
+ "E743", # Ball.I is the moment of inertia, not an ambiguous name
+ "SIM108", # explicit if/else reads better than ternaries in the physics code
+ "B905", # zip(strict=) adds noise to fixed-length numpy pairings
+ "RUF002", # docstrings use unicode maths symbols on purpose
+ "RUF003",
+]
+
+[tool.ruff.lint.per-file-ignores]
+"tests/*" = ["E501"]
+
+[tool.mypy]
+# No python_version pin: mypy targets whichever interpreter runs it, so each entry
+# in the CI matrix checks its own version against the numpy stubs built for it.
+# Ruff's py310 target is what holds the floor for our own syntax.
+files = ["src/cueai"]
+ignore_missing_imports = true
+warn_unused_ignores = true
+warn_redundant_casts = true
diff --git a/requirements.txt b/requirements.txt
index af7e1be..3e2d241 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,17 +1,13 @@
+# Runtime dependencies. For the desktop UI, OpenCV overlays and dev tooling use
+# the extras instead: pip install -e ".[ui,vision,dev]"
numpy>=1.26
pandas>=2.1
-scipy>=1.11
scikit-learn>=1.4
+joblib>=1.3
torch>=2.1
onnx>=1.15
-onnxscript>=0.1.0
onnxruntime>=1.16
-opencv-python-headless>=4.8
-PyQt6>=6.6
fastapi>=0.109
uvicorn[standard]>=0.27
pydantic>=2.5
-httpx>=0.26
-matplotlib>=3.8
tqdm>=4.66
-pytest>=7.4
diff --git a/scripts/benchmark.py b/scripts/benchmark.py
new file mode 100644
index 0000000..cf6d4e9
--- /dev/null
+++ b/scripts/benchmark.py
@@ -0,0 +1,275 @@
+#!/usr/bin/env python3
+"""
+Measure prediction latency and accuracy, and write docs/BENCHMARKS.md.
+
+Everything the README claims about speed comes from this script. Run it after
+``cueai-train`` so the accuracy section can read models/metrics.json:
+
+ python scripts/benchmark.py
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import platform
+import statistics
+import sys
+import time
+from pathlib import Path
+
+import numpy as np
+
+ROOT = Path(__file__).resolve().parents[1]
+sys.path.insert(0, str(ROOT / "src"))
+
+from cueai.ml.infer import TrajectoryPredictor # noqa: E402
+from cueai.physics.analytic import predict_endpoint # noqa: E402
+from cueai.physics.constants import ShotParams, TableParams # noqa: E402
+from cueai.physics.simulator import Simulator # noqa: E402
+
+
+def time_calls(fn, repeats: int) -> dict[str, float]:
+ """Wall-clock statistics for a single call, in milliseconds."""
+ samples = []
+ for _ in range(repeats):
+ start = time.perf_counter()
+ fn()
+ samples.append((time.perf_counter() - start) * 1000)
+ return {
+ "mean_ms": statistics.mean(samples),
+ "median_ms": statistics.median(samples),
+ "p95_ms": float(np.percentile(samples, 95)),
+ "n": repeats,
+ }
+
+
+def measure_latency(repeats: int) -> dict[str, dict[str, float]]:
+ table = TableParams()
+ shot = ShotParams(speed=2.5, angle=0.12, english_x=0.2, english_y=-0.1)
+ cue = (0.6, 0.635)
+ results: dict[str, dict[str, float]] = {}
+
+ rack_sim = Simulator(table=table, dt=0.001, max_time=15.0)
+ results["simulator_full_rack"] = time_calls(
+ lambda: rack_sim.simulate_shot(shot, full_rack=True), max(repeats // 200, 3)
+ )
+
+ two_ball = Simulator(table=table, dt=0.002, max_time=8.0)
+ results["simulator_two_ball"] = time_calls(
+ lambda: two_ball.simulate_shot(shot, cue_pos=cue, obj_pos=(1.4, 0.7)),
+ max(repeats // 50, 5),
+ )
+
+ results["closed_form"] = time_calls(
+ lambda: predict_endpoint(shot, cue, table), repeats
+ )
+
+ gbm_path = ROOT / "models" / "gbm_baseline.joblib"
+ if gbm_path.exists():
+ import joblib
+
+ from cueai.ml.features import build_features
+
+ bundle = joblib.load(gbm_path)
+ row = build_features(shot, cue, (1.4, 0.7), table)[None, :]
+ results["gradient_boosting"] = time_calls(
+ lambda: bundle["model"].predict(bundle["scaler"].transform(row)), repeats
+ )
+
+ predictor = TrajectoryPredictor(model_dir=ROOT / "models")
+ if predictor.ready:
+ results["surrogate_" + predictor.backend] = time_calls(
+ lambda: predictor.predict_fast(shot, cue, obj_pos=(1.4, 0.7), table=table),
+ repeats,
+ )
+ batch = np.repeat(
+ predictor.feature_vector(shot, cue, (1.4, 0.7), table)[None, :], 1024, axis=0
+ )
+ batched = time_calls(lambda: predictor.residual_batch(batch), max(repeats // 20, 5))
+ results["cuenet_batch1024"] = {
+ key: value / 1024 if key.endswith("_ms") else value
+ for key, value in batched.items()
+ }
+
+ reference = results["simulator_full_rack"]["mean_ms"]
+ for name, entry in results.items():
+ if name != "simulator_full_rack":
+ entry["speedup_vs_full_rack"] = reference / entry["mean_ms"]
+ return results
+
+
+def format_duration(milliseconds: float) -> str:
+ """Human-readable across the five orders of magnitude this table spans."""
+ if milliseconds >= 1000:
+ return f"{milliseconds / 1000:,.1f} s"
+ if milliseconds >= 1:
+ return f"{milliseconds:.1f} ms"
+ if milliseconds >= 0.01:
+ return f"{milliseconds:.3f} ms"
+ return f"{milliseconds * 1000:.1f} us"
+
+
+def render_markdown(latency: dict, metrics: dict | None) -> str:
+ labels = {
+ "simulator_full_rack": "Numerical simulator, 16-ball rack",
+ "simulator_two_ball": "Numerical simulator, cue + object ball",
+ "closed_form": "Closed-form solver (no ML)",
+ "gradient_boosting": "Gradient boosting on the same features",
+ "surrogate_torch": "Closed form + CueNet residual (PyTorch)",
+ "surrogate_onnx": "Closed form + CueNet residual (ONNX Runtime)",
+ "cuenet_batch1024": "CueNet forward pass only, batch of 1024 (per shot)",
+ }
+ lines = [
+ "# Benchmarks",
+ "",
+ "Regenerate with `python scripts/benchmark.py`. Single CPU core, no GPU.",
+ "",
+ f"- Python {platform.python_version()} on {platform.platform()}",
+ f"- Processor: {platform.processor() or 'unknown'}",
+ "",
+ "## Latency per shot",
+ "",
+ "| Method | Mean | Median | p95 | Speedup vs full rack |",
+ "|---|---:|---:|---:|---:|",
+ ]
+ for key, entry in latency.items():
+ # The batched row times the network alone, so it is not a like-for-like
+ # end-to-end prediction and does not get a speedup figure.
+ speedup = entry.get("speedup_vs_full_rack") if "batch" not in key else None
+ lines.append(
+ f"| {labels.get(key, key)} | {format_duration(entry['mean_ms'])} | "
+ f"{format_duration(entry['median_ms'])} | {format_duration(entry['p95_ms'])} | "
+ f"{f'{speedup:,.0f}x' if speedup else '—'} |"
+ )
+
+ if metrics:
+ models = metrics["models"]
+ lines += [
+ "",
+ "## Accuracy on held-out shots",
+ "",
+ f"{metrics['n_test']:,} test shots from a {metrics['n_samples']:,} shot dataset. "
+ "Error is the distance between predicted and simulated resting position, "
+ "averaged over the cue ball and the object ball.",
+ "",
+ "| Model | Mean error | p95 | Cue ball | Object ball | R² |",
+ "|---|---:|---:|---:|---:|---:|",
+ ]
+ names = {
+ "analytic": "Closed form (no fitting)",
+ "gbm": "Gradient boosting on raw features",
+ "cuenet": "Closed form + CueNet residual",
+ }
+ for key, label in names.items():
+ m = models[key]
+ lines.append(
+ f"| {label} | {m['euclidean_mm']:.0f} mm | {m['p95_mm']:.0f} mm | "
+ f"{m['cue_mm']:.0f} mm | {m['obj_mm']:.0f} mm | {m['r2']:.3f} |"
+ )
+ lines += [
+ "",
+ "## Where prediction stops working",
+ "",
+ "Resting position is a smooth function of the shot until the ball starts "
+ "ricocheting between cushions. Mean error in mm by cushion contacts:",
+ "",
+ "| Cushion contacts | Shots | Closed form | Gradient boosting | CueNet residual |",
+ "|---|---:|---:|---:|---:|",
+ ]
+ for row in metrics.get("by_cushion_contacts", []):
+ lines.append(
+ f"| {row['cushion_contacts']} | {row['n']} | {row['analytic']:.0f} | "
+ f"{row['gbm']:.0f} | {row['cuenet']:.0f} |"
+ )
+
+ if metrics and metrics.get("by_object_ball_contact"):
+ lines += [
+ "",
+ "## Reading the object-ball number honestly",
+ "",
+ "Only about a third of sampled shots actually hit the object ball. For the "
+ "rest it stays where it started, which the closed-form baseline predicts "
+ "exactly, so those rows contribute almost no object-ball error to any model. "
+ "Split out, in mm:",
+ "",
+ "| Ball-ball contact | Shots | Share | Closed form cue / object | "
+ "Gradient boosting cue / object | CueNet cue / object |",
+ "|---|---:|---:|---:|---:|---:|",
+ ]
+ for row in metrics["by_object_ball_contact"]:
+ cells = " | ".join(
+ f"{row[f'{model}_cue_mm']:.0f} / {row[f'{model}_obj_mm']:.0f}"
+ for model in ("analytic", "gbm", "cuenet")
+ )
+ lines.append(
+ f"| {row['object_ball_contact']} | {row['n']} | "
+ f"{row['share_pct']:.1f}% | {cells} |"
+ )
+
+ if metrics and metrics.get("feature_ablation"):
+ ablation = metrics["feature_ablation"]
+ lines += [
+ "",
+ "## Feature ablation",
+ "",
+ "Same architecture, epochs, seed and split; the ablated model sees only raw "
+ "shot parameters, with none of the closed-form solver's conclusions. This is "
+ "the measurement behind the claim that the features carry this result:",
+ "",
+ "| CueNet inputs | All shots | Direct shots (no cushion) |",
+ "|---|---:|---:|",
+ f"| Raw shot parameters only | {ablation['shot_features_only_mm']:.0f} mm | "
+ f"{ablation['shot_features_only_direct_mm']:.0f} mm |",
+ f"| Plus closed-form output | {ablation['with_closed_form_features_mm']:.0f} mm | "
+ f"{ablation['with_closed_form_features_direct_mm']:.0f} mm |",
+ f"| _closed form alone, for reference_ | "
+ f"{metrics['models']['analytic']['euclidean_mm']:.0f} mm | "
+ f"{ablation['analytic_direct_mm']:.0f} mm |",
+ ]
+
+ if metrics and metrics.get("risk_coverage"):
+ lines += [
+ "",
+ "## Knowing which predictions to trust",
+ "",
+ "The cushion-contact breakdown above slices by what the simulator did, which "
+ "is only knowable after paying for the simulation. This table slices by the "
+ "cushion count the closed-form solver *expects*, which is already computed as "
+ "part of making the prediction and therefore costs nothing extra. Answering "
+ "only the shots below a threshold trades coverage for accuracy:",
+ "",
+ "| Expected cushions | Shots answered | Coverage | Closed form | "
+ "Gradient boosting | CueNet residual |",
+ "|---|---:|---:|---:|---:|---:|",
+ ]
+ for row in metrics["risk_coverage"]:
+ gate = row["expected_cushions_at_most"]
+ label = "no gate" if gate == "all" else f"≤ {gate}"
+ lines.append(
+ f"| {label} | {row['n']} | {row['coverage_pct']:.1f}% | "
+ f"{row['analytic']:.0f} mm | {row['gbm']:.0f} mm | {row['cuenet']:.0f} mm |"
+ )
+ return "\n".join(lines) + "\n"
+
+
+def main() -> None:
+ parser = argparse.ArgumentParser(description=__doc__)
+ parser.add_argument("--repeats", type=int, default=500)
+ parser.add_argument("--out", type=str, default="docs/BENCHMARKS.md")
+ args = parser.parse_args()
+
+ latency = measure_latency(args.repeats)
+ metrics_path = ROOT / "models" / "metrics.json"
+ metrics = json.loads(metrics_path.read_text()) if metrics_path.exists() else None
+
+ (ROOT / "models" / "latency.json").write_text(json.dumps(latency, indent=2) + "\n")
+ out_path = ROOT / args.out
+ out_path.parent.mkdir(parents=True, exist_ok=True)
+ out_path.write_text(render_markdown(latency, metrics))
+ print(json.dumps(latency, indent=2))
+ print(f"wrote {out_path.relative_to(ROOT)}")
+
+
+if __name__ == "__main__":
+ main()
diff --git a/scripts/export_parity_cases.py b/scripts/export_parity_cases.py
new file mode 100644
index 0000000..6a738f6
--- /dev/null
+++ b/scripts/export_parity_cases.py
@@ -0,0 +1,338 @@
+#!/usr/bin/env python3
+"""
+Export reference shots so the browser physics can be checked against Python.
+
+The Python package under `src/cueai/physics/` is the definition of correct: it
+is what `tests/test_validation.py` pins to closed-form mechanics. The browser
+runs a hand port of it, and a port is only worth anything if someone measures
+the difference. This writes the initial conditions and the reference outcome
+for a spread of shots; `web/test/parity.mjs` replays them in Node.
+
+Each case also carries a chaos yardstick. Break shots are Lyapunov unstable, so
+two implementations that differ in the last bit of a float will not agree on
+where sixteen balls come to rest, and demanding that they do would be a
+misunderstanding rather than a standard. The yardstick is how far the reference
+moves when its own initial condition is nudged by a picometre, which bounds how
+much agreement is available to ask for.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import time
+from pathlib import Path
+
+import numpy as np
+
+from cueai.physics.ball import Ball
+from cueai.physics.constants import BallParams, ShotParams, TableParams
+from cueai.physics.rack import identity_for, make_full_rack
+from cueai.physics.simulator import Simulator
+
+PERTURBATION_M = 1e-12
+
+
+def make_ball(number: int, pos: tuple[float, float], params: BallParams) -> Ball:
+ return Ball(
+ id=number,
+ number=number,
+ pos=np.array(pos, dtype=np.float64),
+ vel=np.zeros(2),
+ omega=np.zeros(3),
+ params=params,
+ identity=identity_for(number),
+ )
+
+
+def scatter_layout(
+ rng: np.random.Generator, table: TableParams, params: BallParams, n_objects: int
+) -> list[Ball]:
+ """Place a cue ball and a few object balls with no initial overlap."""
+ R = params.radius
+ margin = 3 * R
+ placed: list[tuple[float, float]] = []
+ while len(placed) < n_objects + 1:
+ candidate = (
+ float(rng.uniform(margin, table.length - margin)),
+ float(rng.uniform(margin, table.width - margin)),
+ )
+ if any(np.hypot(candidate[0] - x, candidate[1] - y) < 4 * R for x, y in placed):
+ continue
+ # Keep them off the pocket mouths so the layout itself is not a pot.
+ if any(np.hypot(candidate[0] - px, candidate[1] - py) < 3 * R for px, py in table.pockets):
+ continue
+ placed.append(candidate)
+
+ numbers = [0, *rng.choice(np.arange(1, 16), size=n_objects, replace=False)]
+ return [make_ball(int(n), pos, params) for n, pos in zip(numbers, placed)]
+
+
+def serialize_balls(balls: list[Ball]) -> list[dict]:
+ return [
+ {"number": int(b.number), "x": float(b.pos[0]), "y": float(b.pos[1])} for b in balls
+ ]
+
+
+def outcome(sim: Simulator, balls: list[Ball], shot: ShotParams) -> dict:
+ started = time.perf_counter()
+ result = sim.simulate_shot(shot, balls=balls)
+ elapsed = time.perf_counter() - started
+ return {
+ "seconds": elapsed,
+ "table_time": float(result.times[-1]) if len(result.times) else 0.0,
+ "resting": {
+ str(int(b.number)): [float(result.endpoints[b.id][0]), float(result.endpoints[b.id][1])]
+ for b in balls
+ if not result.pocketed[b.id]
+ },
+ "pocketed": sorted(int(b.number) for b in balls if result.pocketed[b.id]),
+ "collisions": int(result.collision_events),
+ "cushions": int(result.cushion_events),
+ }
+
+
+def chaos_yardstick(sim: Simulator, balls: list[Ball], shot: ShotParams, reference: dict) -> float:
+ """
+ How far the reference moves when its own input is nudged by a picometre.
+
+ A porting difference smaller than this is indistinguishable from the
+ sensitivity the system already has to the last bit of its inputs.
+ """
+ nudged = [b.copy() for b in balls]
+ nudged[0].pos[0] += PERTURBATION_M
+ other = outcome(sim, nudged, shot)
+ worst = 0.0
+ for number, pos in reference["resting"].items():
+ if number not in other["resting"]:
+ return float("inf") # a ball changed pocketed status: fully chaotic
+ elsewhere = other["resting"][number]
+ moved = float(np.hypot(pos[0] - elsewhere[0], pos[1] - elsewhere[1]))
+ worst = max(worst, moved)
+ if reference["pocketed"] != other["pocketed"]:
+ return float("inf")
+ return worst
+
+
+def build_cases(seed: int) -> list[dict]:
+ rng = np.random.default_rng(seed)
+ table = TableParams()
+ params = BallParams()
+ sim = Simulator(table=table, ball_params=params)
+ cases: list[dict] = []
+
+ # Hand-picked shots that exercise one mechanism each, so a parity failure
+ # points at the part of the model that was ported wrong.
+ named: list[tuple[str, list[Ball], ShotParams]] = [
+ (
+ "stun-into-rail",
+ [make_ball(0, (0.6, 0.635), params)],
+ ShotParams(speed=3.0, angle=0.0),
+ ),
+ (
+ "draw",
+ [make_ball(0, (0.6, 0.635), params), make_ball(1, (1.6, 0.635), params)],
+ ShotParams(speed=3.5, angle=0.0, english_y=-0.45),
+ ),
+ (
+ "follow",
+ [make_ball(0, (0.6, 0.635), params), make_ball(1, (1.6, 0.635), params)],
+ ShotParams(speed=3.5, angle=0.0, english_y=0.45),
+ ),
+ (
+ "right-english-off-two-rails",
+ [make_ball(0, (0.6, 0.4), params)],
+ ShotParams(speed=4.5, angle=0.7, english_x=0.45),
+ ),
+ (
+ "thin-cut",
+ [make_ball(0, (0.6, 0.5), params), make_ball(1, (1.5, 0.7), params)],
+ ShotParams(speed=4.0, angle=0.2),
+ ),
+ (
+ "corner-pocket",
+ [make_ball(0, (1.2, 0.635), params), make_ball(1, (2.0, 0.9), params)],
+ ShotParams(speed=3.0, angle=float(np.arctan2(0.9 - 0.635, 2.0 - 1.2))),
+ ),
+ (
+ "three-ball-cluster",
+ [
+ make_ball(0, (0.5, 0.635), params),
+ make_ball(1, (1.5, 0.635), params),
+ make_ball(2, (1.5 + 2 * params.radius + 1e-4, 0.66), params),
+ make_ball(3, (1.5 + 2 * params.radius + 1e-4, 0.60), params),
+ ],
+ ShotParams(speed=5.0, angle=0.0),
+ ),
+ (
+ "soft-roll",
+ [make_ball(0, (0.5, 0.635), params)],
+ ShotParams(speed=0.8, angle=0.35),
+ ),
+ (
+ "masse-lite",
+ [make_ball(0, (0.9, 0.5), params)],
+ ShotParams(speed=2.5, angle=1.2, english_x=-0.4, cue_elevation=0.12),
+ ),
+ ]
+
+ for name, balls, shot in named:
+ ref = outcome(sim, balls, shot)
+ cases.append(
+ {
+ "name": name,
+ "balls": serialize_balls(balls),
+ "shot": {
+ "speed": shot.speed,
+ "angle": shot.angle,
+ "english_x": shot.english_x,
+ "english_y": shot.english_y,
+ "cue_elevation": shot.cue_elevation,
+ },
+ "reference": ref,
+ "chaos_yardstick_m": chaos_yardstick(sim, balls, shot, ref),
+ }
+ )
+
+ # Randomised layouts, to catch anything the hand-picked shots miss.
+ for i in range(24):
+ balls = scatter_layout(rng, table, params, n_objects=int(rng.integers(1, 5)))
+ shot = ShotParams(
+ speed=float(rng.uniform(1.0, 6.0)),
+ angle=float(rng.uniform(-np.pi, np.pi)),
+ english_x=float(rng.uniform(-0.45, 0.45)),
+ english_y=float(rng.uniform(-0.45, 0.45)),
+ )
+ ref = outcome(sim, balls, shot)
+ cases.append(
+ {
+ "name": f"random-{i:02d}",
+ "balls": serialize_balls(balls),
+ "shot": {
+ "speed": shot.speed,
+ "angle": shot.angle,
+ "english_x": shot.english_x,
+ "english_y": shot.english_y,
+ "cue_elevation": 0.0,
+ },
+ "reference": ref,
+ "chaos_yardstick_m": chaos_yardstick(sim, balls, shot, ref),
+ }
+ )
+
+ # The break: sixteen balls, the case the game actually opens with.
+ for speed in (6.0, 8.0):
+ balls = make_full_rack(table=table, ball_params=params, seed=7)
+ shot = ShotParams(speed=speed, angle=0.0)
+ ref = outcome(sim, balls, shot)
+ cases.append(
+ {
+ "name": f"break-{speed:.0f}ms",
+ "balls": serialize_balls(balls),
+ "shot": {
+ "speed": shot.speed,
+ "angle": shot.angle,
+ "english_x": 0.0,
+ "english_y": 0.0,
+ "cue_elevation": 0.0,
+ },
+ "reference": ref,
+ "chaos_yardstick_m": chaos_yardstick(sim, balls, shot, ref),
+ }
+ )
+
+ return cases
+
+
+def physics_only(case: dict) -> dict:
+ """The reproducible part of a case, with wall-clock timings dropped."""
+ return {
+ "name": case["name"],
+ "balls": case["balls"],
+ "shot": case["shot"],
+ "pocketed": case["reference"]["pocketed"],
+ "resting": case["reference"]["resting"],
+ }
+
+
+def check_against(path: Path, cases: list[dict]) -> int:
+ """
+ Verify the committed cases still describe what the simulator does.
+
+ Timings vary run to run, so comparing the files byte for byte would fail
+ for no reason; only the physics is compared.
+ """
+ if not path.exists():
+ print(f"{path} does not exist; run without --check to create it")
+ return 1
+
+ committed = {c["name"]: physics_only(c) for c in json.loads(path.read_text())["cases"]}
+ fresh = {c["name"]: physics_only(c) for c in cases}
+
+ if committed.keys() != fresh.keys():
+ print(f"case list changed: {sorted(set(fresh) ^ set(committed))}")
+ return 1
+
+ worst = 0.0
+ worst_name = ""
+ for name, current in fresh.items():
+ old = committed[name]
+ if old["pocketed"] != current["pocketed"]:
+ print(f"{name}: pocketed {old['pocketed']} but now pockets {current['pocketed']}")
+ return 1
+ for number, pos in current["resting"].items():
+ before = old["resting"][number]
+ moved = float(np.hypot(pos[0] - before[0], pos[1] - before[1]))
+ if moved > worst:
+ worst, worst_name = moved, f"{name} ball {number}"
+
+ if worst > 1e-6:
+ print(f"committed cases are stale: {worst * 1000:.4f} mm drift on {worst_name}")
+ print("run `make parity` and commit web/test/parity_cases.json")
+ return 1
+ print(f"committed cases match the simulator (worst drift {worst * 1000:.2e} mm)")
+ return 0
+
+
+def main(argv: list[str] | None = None) -> None:
+ parser = argparse.ArgumentParser(description=__doc__)
+ parser.add_argument("--out", type=Path, default=Path("web/test/parity_cases.json"))
+ parser.add_argument("--seed", type=int, default=11)
+ parser.add_argument(
+ "--check",
+ action="store_true",
+ help="compare against the committed cases instead of rewriting them",
+ )
+ args = parser.parse_args(argv)
+
+ cases = build_cases(args.seed)
+
+ if args.check:
+ raise SystemExit(check_against(args.out, cases))
+
+ args.out.parent.mkdir(parents=True, exist_ok=True)
+ payload = {
+ "dt": 0.001,
+ "max_time": 15.0,
+ "perturbation_m": PERTURBATION_M,
+ "note": (
+ "Generated by scripts/export_parity_cases.py from the Python reference "
+ "simulator. chaos_yardstick_m is how far that reference moves under a "
+ "1e-12 m nudge to the cue ball, which bounds the agreement any second "
+ "implementation can be asked for."
+ ),
+ "cases": cases,
+ }
+ args.out.write_text(json.dumps(payload, indent=1) + "\n")
+
+ finite = [c["chaos_yardstick_m"] for c in cases if np.isfinite(c["chaos_yardstick_m"])]
+ total = sum(c["reference"]["seconds"] for c in cases)
+ table_time = sum(c["reference"]["table_time"] for c in cases)
+ print(f"wrote {len(cases)} cases to {args.out}")
+ print(f" {len(finite)} deterministic enough to compare tightly")
+ print(f" {len(cases) - len(finite)} chaotic (a picometre changes which balls drop)")
+ print(f" reference spent {total:.2f} s simulating {table_time:.1f} s of table time")
+
+
+if __name__ == "__main__":
+ main()
diff --git a/scripts/make_figures.py b/scripts/make_figures.py
new file mode 100644
index 0000000..9f2edcb
--- /dev/null
+++ b/scripts/make_figures.py
@@ -0,0 +1,331 @@
+#!/usr/bin/env python3
+"""
+Render the figures used in the README, headlessly.
+
+ python scripts/make_figures.py
+
+Writes docs/assets/*.png. The validation and break figures need nothing but the
+physics; the accuracy and latency figures read models/metrics.json and
+models/latency.json, so run training and scripts/benchmark.py first.
+"""
+
+from __future__ import annotations
+
+import json
+import sys
+from pathlib import Path
+
+import matplotlib
+
+matplotlib.use("Agg")
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+ROOT = Path(__file__).resolve().parents[1]
+sys.path.insert(0, str(ROOT / "src"))
+
+from cueai.physics import analytic # noqa: E402
+from cueai.physics.ball import Ball, MotionState, integrate_ball # noqa: E402
+from cueai.physics.constants import BallParams, ShotParams, TableParams # noqa: E402
+from cueai.physics.simulator import Simulator # noqa: E402
+
+ASSETS = ROOT / "docs" / "assets"
+CLOTH = "#12764a"
+INK = "#1b1f23"
+ACCENT = "#e06c1f"
+BLUE = "#2f6fb0"
+
+plt.rcParams.update(
+ {
+ "figure.dpi": 130,
+ "savefig.bbox": "tight",
+ "font.size": 9,
+ "axes.spines.top": False,
+ "axes.spines.right": False,
+ "axes.grid": True,
+ "grid.alpha": 0.25,
+ "grid.linestyle": ":",
+ }
+)
+
+
+def _roll_to_transition(v0: float, table: TableParams) -> tuple[float, float]:
+ """Integrate a centre-ball hit up to the sliding/rolling transition."""
+ params = BallParams()
+ ball = Ball(
+ id=0,
+ pos=np.array([0.0, 0.635]),
+ vel=np.array([v0, 0.0]),
+ omega=np.zeros(3),
+ params=params,
+ )
+ for _ in range(400_000):
+ ball = integrate_ball(ball, table, 1e-4)
+ if ball.motion_state(table) is MotionState.ROLLING:
+ return ball.speed(), float(ball.pos[0])
+ raise RuntimeError("ball never reached the rolling phase")
+
+
+def figure_validation() -> Path:
+ """Simulator versus closed form for the sliding phase."""
+ table = TableParams(friction_noise_amp=0.0)
+ speeds = np.linspace(0.6, 4.0, 12)
+ measured = np.array([_roll_to_transition(float(v), table) for v in speeds])
+
+ fig, axes = plt.subplots(1, 2, figsize=(8.2, 3.2))
+
+ axes[0].plot(speeds, analytic.ROLLING_SPEED_RATIO * speeds, "-", color=INK, lw=1.4,
+ label=r"closed form $\frac{5}{7}v_0$")
+ axes[0].plot(speeds, measured[:, 0], "o", color=ACCENT, ms=4.5, label="simulator")
+ axes[0].set_xlabel("launch speed $v_0$ (m/s)")
+ axes[0].set_ylabel("speed at rolling (m/s)")
+ axes[0].set_title("Sliding ends at 5/7 of launch speed")
+ axes[0].legend(frameon=False)
+
+ reference = [analytic.slide_distance(float(v), table.mu_slide) for v in speeds]
+ axes[1].plot(speeds, reference, "-", color=INK, lw=1.4,
+ label=r"closed form $\frac{12v_0^2}{49\mu_s g}$")
+ axes[1].plot(speeds, measured[:, 1], "o", color=BLUE, ms=4.5, label="simulator")
+ axes[1].set_xlabel("launch speed $v_0$ (m/s)")
+ axes[1].set_ylabel("sliding distance (m)")
+ axes[1].set_title("Sliding distance matches theory")
+ axes[1].legend(frameon=False)
+
+ worst = float(
+ np.max(np.abs(measured[:, 1] - reference) / np.maximum(reference, 1e-9)) * 100
+ )
+ fig.suptitle(
+ f"Physics validation: worst-case deviation from closed form {worst:.1f}%",
+ fontsize=10,
+ )
+ return _save(fig, "validation.png")
+
+
+def _draw_table(ax, table: TableParams) -> None:
+ ax.add_patch(
+ plt.Rectangle((0, 0), table.length, table.width, facecolor=CLOTH, edgecolor="#5c3a1e", lw=6)
+ )
+ for x, y in table.pockets:
+ ax.add_patch(plt.Circle((x, y), table.pocket_radius, color="#0b0b0b", zorder=3))
+ ax.set_xlim(-0.08, table.length + 0.08)
+ ax.set_ylim(-0.08, table.width + 0.08)
+ ax.set_aspect("equal")
+ ax.axis("off")
+
+
+def figure_break() -> Path:
+ """A full 16-ball break, which is what the simulator is actually solving."""
+ table = TableParams(friction_noise_amp=0.02)
+ sim = Simulator(table=table, dt=0.001, max_time=15.0, collision_passes=20)
+ shot = ShotParams(speed=6.5, angle=0.008, english_y=0.25)
+ result = sim.simulate_shot(shot, full_rack=True, seed=7)
+
+ fig, ax = plt.subplots(figsize=(8.6, 4.6))
+ _draw_table(ax, table)
+ for ball_id, path in result.trajectories.items():
+ visible = path[np.all(path >= 0, axis=1)]
+ if len(visible) < 2:
+ continue
+ meta = result.ball_meta.get(ball_id, {})
+ colour = np.array(meta.get("color", [220, 220, 220])) / 255
+ ax.plot(visible[:, 0], visible[:, 1], "-", color=colour, lw=1.6, alpha=0.9, zorder=4)
+ ax.plot(visible[-1, 0], visible[-1, 1], "o", color=colour, ms=6,
+ markeredgecolor="#111", markeredgewidth=0.6, zorder=5)
+
+ potted = sum(1 for was_potted in result.pocketed.values() if was_potted)
+ ax.set_title(
+ f"16-ball break: {result.collision_events} ball-ball contacts, "
+ f"{result.cushion_events} cushion contacts, {potted} potted",
+ fontsize=10,
+ )
+ return _save(fig, "break_shot.png")
+
+
+def figure_spin() -> Path:
+ """Draw, stun and follow from the same stroke speed."""
+ table = TableParams(friction_noise_amp=0.0)
+ sim = Simulator(table=table, dt=5e-4, max_time=12.0, collision_passes=10)
+ fig, ax = plt.subplots(figsize=(8.6, 4.6))
+ _draw_table(ax, table)
+
+ # One lane per stroke so the three cue-ball paths do not overlap.
+ # The gap is kept short so the centre-ball shot is still sliding at contact,
+ # which is what makes the three strokes visibly different.
+ start_x = 0.95
+ lanes = [
+ (0.95, -0.45, "backspin (tip 0.45R below centre) -> draw", ACCENT),
+ (0.635, 0.0, "centre ball -> stun", "#f4f4f4"),
+ (0.32, 0.45, "topspin (tip 0.45R above centre) -> follow", BLUE),
+ ]
+ for lane_y, english_y, label, colour in lanes:
+ result = sim.simulate_shot(
+ ShotParams(speed=2.5, angle=0.0, english_y=english_y),
+ cue_pos=(start_x, lane_y),
+ obj_pos=(1.2, lane_y),
+ full_rack=False,
+ )
+ cue = result.trajectories[0]
+ ax.plot(cue[:, 0], cue[:, 1], "-", color=colour, lw=2.2, zorder=4)
+ ax.plot(start_x, lane_y, "o", color=colour, ms=9, markeredgecolor="#111", zorder=5)
+ ax.plot(cue[-1, 0], cue[-1, 1], "*", color=colour, ms=15, markeredgecolor="#111",
+ markeredgewidth=0.5, zorder=6)
+ ax.plot(1.2, lane_y, "o", color="#e8e2b0", ms=9, markeredgecolor="#111", zorder=5)
+ ax.text(0.10, lane_y + 0.09, label, color=colour, fontsize=8.5, zorder=7)
+ ax.text(float(cue[-1, 0]), lane_y - 0.13, f"{float(cue[-1, 0]):.2f} m",
+ color=colour, fontsize=8, ha="center", zorder=7)
+
+ ax.annotate(
+ "contact",
+ xy=(1.2, 0.05),
+ xytext=(1.2, -0.02),
+ color="#f4f4f4",
+ fontsize=8,
+ ha="center",
+ )
+ ax.axvline(1.2, color="#f4f4f4", alpha=0.25, lw=0.8, ls=":", zorder=2)
+ ax.set_title(
+ "Same 2.5 m/s stroke, three tip heights: the cue ball draws back, stops short, "
+ "or follows through\n(circles = start, stars = resting position)",
+ fontsize=9.5,
+ )
+ return _save(fig, "spin_control.png")
+
+
+def figure_accuracy(metrics: dict) -> Path:
+ """Error by cushion contacts, which is where the approach runs out."""
+ rows = metrics["by_cushion_contacts"]
+ labels = [row["cushion_contacts"] for row in rows]
+ series = [
+ ("closed form", "analytic", "#9aa4ad"),
+ ("gradient boosting", "gbm", BLUE),
+ ("closed form + CueNet", "cuenet", ACCENT),
+ ]
+ positions = np.arange(len(labels))
+ width = 0.26
+
+ fig, ax = plt.subplots(figsize=(7.2, 3.4))
+ for offset, (label, key, colour) in enumerate(series):
+ values = [row[key] for row in rows]
+ ax.bar(positions + (offset - 1) * width, values, width, label=label, color=colour)
+ ax.set_xticks(positions)
+ ax.set_xticklabels([f"{label}\n(n={row['n']})" for label, row in zip(labels, rows)])
+ ax.set_xlabel("cushion contacts before coming to rest")
+ ax.set_ylabel("mean endpoint error (mm)")
+ ax.set_title(
+ "Prediction error grows with every cushion contact, until the outcome is chaotic",
+ fontsize=10,
+ )
+ ax.legend(frameon=False)
+ return _save(fig, "accuracy.png")
+
+
+def figure_reliability(metrics: dict) -> Path:
+ """Error against coverage, when the closed-form cushion count is used as a gate."""
+ rows = metrics["risk_coverage"]
+ coverage = [row["coverage_pct"] for row in rows]
+ series = [
+ ("closed form", "analytic", "#9aa4ad", "o"),
+ ("gradient boosting", "gbm", BLUE, "s"),
+ ("closed form + CueNet", "cuenet", ACCENT, "D"),
+ ]
+
+ fig, ax = plt.subplots(figsize=(7.2, 3.4))
+ for label, key, colour, marker in series:
+ ax.plot(coverage, [row[key] for row in rows], marker + "-", color=colour,
+ lw=1.6, ms=5, label=label)
+ for row in rows:
+ gate = row["expected_cushions_at_most"]
+ ax.annotate(
+ "no gate" if gate == "all" else f"≤{gate}",
+ xy=(row["coverage_pct"], row["cuenet"]),
+ xytext=(0, -13),
+ textcoords="offset points",
+ ha="center",
+ fontsize=7.5,
+ color=ACCENT,
+ )
+ ax.set_xlabel("shots answered by the fast path (%)")
+ ax.set_ylabel("mean endpoint error (mm)")
+ ax.set_title(
+ "Choosing how much to answer: the closed-form cushion count is a free gate",
+ fontsize=10,
+ )
+ ax.legend(frameon=False, loc="upper left")
+ return _save(fig, "reliability.png")
+
+
+def figure_latency(latency: dict) -> Path:
+ """Cost per shot, log scale, because the range spans five orders of magnitude."""
+ labels = {
+ "simulator_full_rack": "simulator, 16 balls",
+ "simulator_two_ball": "simulator, 2 balls",
+ "closed_form": "closed form",
+ "surrogate_onnx": "closed form + CueNet (ONNX)",
+ "surrogate_torch": "closed form + CueNet (PyTorch)",
+ "cuenet_batch1024": "CueNet forward pass, batched",
+ }
+ entries = [
+ (labels[key], value["mean_ms"]) for key, value in latency.items() if key in labels
+ ]
+ entries.sort(key=lambda item: item[1])
+
+ def label_time(milliseconds: float) -> str:
+ if milliseconds >= 1000:
+ return f"{milliseconds / 1000:,.1f} s"
+ if milliseconds >= 1:
+ return f"{milliseconds:.0f} ms"
+ if milliseconds >= 0.01:
+ return f"{milliseconds:.2f} ms"
+ return f"{milliseconds * 1000:.1f} us"
+
+ fig, ax = plt.subplots(figsize=(7.2, 3.0))
+ names = [name for name, _ in entries]
+ values = [value for _, value in entries]
+ colours = [ACCENT if "CueNet" in name or "closed" in name else "#9aa4ad" for name in names]
+ ax.barh(names, values, color=colours)
+ ax.set_xscale("log")
+ ax.set_xlabel("time per shot, log scale")
+ for index, value in enumerate(values):
+ ax.text(value * 1.15, index, label_time(value), va="center", fontsize=8)
+ ax.set_xlim(min(values) * 0.5, max(values) * 6)
+ ax.set_title("Cost of one shot prediction", fontsize=10)
+ ax.grid(axis="y", visible=False)
+ return _save(fig, "latency.png")
+
+
+def _save(fig, name: str) -> Path:
+ ASSETS.mkdir(parents=True, exist_ok=True)
+ path = ASSETS / name
+ fig.savefig(path)
+ plt.close(fig)
+ print(f"wrote {path.relative_to(ROOT)}")
+ return path
+
+
+def main() -> None:
+ figure_validation()
+ figure_spin()
+ figure_break()
+
+ metrics_path = ROOT / "models" / "metrics.json"
+ metrics = json.loads(metrics_path.read_text()) if metrics_path.exists() else {}
+ if metrics.get("by_cushion_contacts"):
+ figure_accuracy(metrics)
+ else:
+ print("skipping accuracy figure: run training to produce models/metrics.json")
+
+ if metrics.get("risk_coverage"):
+ figure_reliability(metrics)
+ else:
+ print("skipping reliability figure: run training to produce models/metrics.json")
+
+ latency_path = ROOT / "models" / "latency.json"
+ if latency_path.exists():
+ figure_latency(json.loads(latency_path.read_text()))
+ else:
+ print("skipping latency figure: models/latency.json missing")
+
+
+if __name__ == "__main__":
+ main()
diff --git a/scripts/site_facts.py b/scripts/site_facts.py
new file mode 100644
index 0000000..82242be
--- /dev/null
+++ b/scripts/site_facts.py
@@ -0,0 +1,195 @@
+"""
+Collect every number the playable page quotes into one measured file.
+
+The explainer on the table page makes specific claims — how many physics tests
+there are, how closely the browser port tracks the reference, what the learned
+surrogate costs and buys. Numbers typed into HTML are true on the day they are
+typed. These are read from the artefacts the test suite, the parity harness and
+the training run leave behind, so the page is either current or visibly stale.
+
+ python scripts/site_facts.py [--check]
+
+``--check`` fails instead of writing, for continuous integration.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import re
+import subprocess
+import sys
+from pathlib import Path
+from typing import Any
+
+ROOT = Path(__file__).resolve().parents[1]
+OUT = ROOT / "web" / "data" / "facts.json"
+PAGE = ROOT / "web" / "index.html"
+TIMES = "\N{MULTIPLICATION SIGN}"
+
+# Wall-clock measurements are a property of the machine that took them, so
+# ``--check`` reports drift in these but does not fail on it. Everything else
+# is a deterministic function of the code and has to match.
+VOLATILE = frozenset(
+ {
+ "parity-speedup",
+ "parity-python-seconds",
+ "parity-browser-seconds",
+ "sim-full-rack",
+ "surrogate-latency",
+ "closed-form-latency",
+ }
+)
+
+
+def _load(path: Path) -> dict[str, Any] | None:
+ try:
+ return json.loads(path.read_text())
+ except (OSError, json.JSONDecodeError):
+ return None
+
+
+def count_physics_tests() -> int | None:
+ """Ask pytest how many validation tests there are, rather than guessing."""
+ try:
+ proc = subprocess.run(
+ [sys.executable, "-m", "pytest", "--collect-only", "-q", "tests/test_validation.py"],
+ cwd=ROOT,
+ capture_output=True,
+ text=True,
+ timeout=180,
+ )
+ except (OSError, subprocess.SubprocessError):
+ return None
+ if proc.returncode != 0:
+ return None
+ for line in reversed(proc.stdout.splitlines()):
+ if "test" in line and "collected" in line:
+ head = line.split()[0]
+ if head.isdigit():
+ return int(head)
+ return None
+
+
+def mm(value: float) -> str:
+ return f"{value:.0f} mm"
+
+
+def build() -> dict[str, str]:
+ facts: dict[str, str] = {}
+
+ tests = count_physics_tests()
+ if tests is not None:
+ facts["physics-tests"] = str(tests)
+
+ parity = _load(ROOT / "web" / "data" / "parity.json")
+ if parity:
+ worst = float(parity["worst_mm"])
+ # Below a micrometre the reader wants the order of magnitude, not the
+ # digits; above it, a plain number reads better than scientific form.
+ mantissa, exponent = f"{worst:.1e}".split("e")
+ facts["parity-worst"] = (
+ f"{mantissa} {TIMES} 10{int(exponent)} mm"
+ if worst < 1e-3
+ else f"{worst:.4f} mm"
+ )
+ facts["parity-cases"] = str(parity["cases"])
+ facts["parity-speedup"] = f"{float(parity['speedup']):.0f}{TIMES}"
+ facts["parity-table-seconds"] = f"{float(parity['table_seconds']):.0f} s"
+ facts["parity-python-seconds"] = f"{float(parity['python_seconds']):.0f} s"
+ facts["parity-browser-seconds"] = f"{float(parity['browser_seconds']):.2f} s"
+
+ metrics = _load(ROOT / "models" / "metrics.json")
+ if metrics:
+ models = metrics["models"]
+ facts["surrogate-error"] = mm(models["cuenet"]["euclidean_mm"])
+ facts["analytic-error"] = mm(models["analytic"]["euclidean_mm"])
+ facts["gbm-error"] = mm(models["gbm"]["euclidean_mm"])
+ facts["samples"] = f"{int(metrics['n_samples']):,}"
+ contact = {row["object_ball_contact"]: row for row in metrics["by_object_ball_contact"]}
+ if "yes" in contact:
+ facts["contact-cuenet"] = mm(contact["yes"]["cuenet_obj_mm"])
+ facts["contact-gbm"] = mm(contact["yes"]["gbm_obj_mm"])
+
+ latency = _load(ROOT / "models" / "latency.json")
+ if latency:
+ full = float(latency["simulator_full_rack"]["mean_ms"])
+ facts["sim-full-rack"] = f"{full / 1000:.1f} s" if full >= 1000 else f"{full:.0f} ms"
+ facts["surrogate-latency"] = f"{float(latency['surrogate_onnx']['mean_ms']):.2f} ms"
+ facts["closed-form-latency"] = f"{float(latency['closed_form']['mean_ms']):.2f} ms"
+
+ return facts
+
+
+def rewrite_page(html: str, facts: dict[str, str]) -> str:
+ """
+ Put the measured values into the page's own markup.
+
+ The loader replaces them at runtime too, but only once the fetch resolves,
+ and never at all over ``file://``. Writing them in as well means the text
+ is never briefly wrong and never wrong at rest. Assumes a ``data-fact``
+ element contains no further element of its own tag, which is true here.
+ """
+
+ def substitute(match: re.Match[str]) -> str:
+ tag, attributes, key = match.group(1), match.group(2), match.group(3)
+ if key not in facts:
+ return match.group(0)
+ return f"<{tag}{attributes}>{facts[key]}{tag}>"
+
+ pattern = re.compile(
+ r"<(\w+)([^>]*\bdata-fact=\"([\w-]+)\"[^>]*)>.*?\1>",
+ re.DOTALL,
+ )
+ return pattern.sub(substitute, html)
+
+
+def main(argv: list[str] | None = None) -> None:
+ parser = argparse.ArgumentParser(description=__doc__)
+ parser.add_argument("--check", action="store_true", help="fail if the file is out of date")
+ args = parser.parse_args(argv)
+
+ facts = build()
+ if not facts:
+ raise SystemExit("no measurements found; run `make all` and `make parity` first")
+
+ payload = f"{json.dumps(facts, indent=2)}\n"
+ html = PAGE.read_text()
+ # Held to the reproducible figures only: timings differ from machine to
+ # machine, so failing on those would fail for a reason nobody can fix.
+ checkable = rewrite_page(html, {k: v for k, v in facts.items() if k not in VOLATILE})
+
+ if args.check:
+ rewritten = checkable
+ committed = _load(OUT) or {}
+ stale = [
+ f" {key}: committed {committed.get(key, '(missing)')!r}, measured {value!r}"
+ for key, value in facts.items()
+ if key not in VOLATILE and committed.get(key) != value
+ ]
+ if stale:
+ print("web/data/facts.json is out of date; run `python scripts/site_facts.py`")
+ print("\n".join(stale))
+ raise SystemExit(1)
+ if rewritten != html:
+ print("web/index.html quotes numbers the measurements no longer support;")
+ print("run `python scripts/site_facts.py`")
+ raise SystemExit(1)
+ drifted = sum(committed.get(k) != v for k, v in facts.items() if k in VOLATILE)
+ print(
+ f"web/data/facts.json is current ({len(facts) - len(VOLATILE)} fixed measurements"
+ + (f", {drifted} timing figure(s) differ on this machine)" if drifted else ")")
+ )
+ return
+
+ OUT.parent.mkdir(parents=True, exist_ok=True)
+ OUT.write_text(payload)
+ rewritten = rewrite_page(html, facts)
+ if rewritten != html:
+ PAGE.write_text(rewritten)
+ print(f"rewrote the figures quoted in {PAGE.relative_to(ROOT)}")
+ print(f"wrote {len(facts)} measurements to {OUT.relative_to(ROOT)}")
+
+
+if __name__ == "__main__":
+ main()
diff --git a/src/cueai/api/main.py b/src/cueai/api/main.py
index 0d07571..999bc1a 100644
--- a/src/cueai/api/main.py
+++ b/src/cueai/api/main.py
@@ -2,20 +2,29 @@
from __future__ import annotations
+import time
from pathlib import Path
+import numpy as np
from fastapi import FastAPI
from fastapi.middleware.cors import CORSMiddleware
from pydantic import BaseModel, Field
+from cueai.ml.dataset import MAX_TIP_OFFSET
from cueai.ml.infer import TrajectoryPredictor
from cueai.physics.constants import ShotParams, TableParams
from cueai.physics.rack import make_full_rack
+API_VERSION = "0.3.0"
+
app = FastAPI(
title="CueAI API",
- description="Physics-informed full-rack billiards simulation",
- version="0.2.0",
+ description=(
+ "Physics-informed billiards prediction. /predict runs the full numerical "
+ "simulation; /predict/fast returns the closed-form plus learned-residual "
+ "estimate in under a millisecond."
+ ),
+ version=API_VERSION,
)
app.add_middleware(
CORSMiddleware,
@@ -31,8 +40,14 @@
class ShotRequest(BaseModel):
speed: float = Field(4.5, ge=0.1, le=12)
angle_deg: float = Field(0.0, description="Launch angle in degrees")
- english_x: float = Field(0.0, ge=-1, le=1)
- english_y: float = Field(0.0, ge=-1, le=1)
+ # Tip offsets are bounded by the miscue limit, which is also where the
+ # training distribution stops: a request outside it would be extrapolation.
+ english_x: float = Field(
+ 0.0, ge=-MAX_TIP_OFFSET, le=MAX_TIP_OFFSET, description="Sidespin tip offset, fraction of R"
+ )
+ english_y: float = Field(
+ 0.0, ge=-MAX_TIP_OFFSET, le=MAX_TIP_OFFSET, description="Top/backspin tip offset"
+ )
cue_elevation_deg: float = 0.0
cue_x: float = 0.635
cue_y: float = 0.635
@@ -48,34 +63,37 @@ class ShotRequest(BaseModel):
class HealthResponse(BaseModel):
status: str
ml_loaded: bool
+ backend: str
version: str
+def _to_physics(req: ShotRequest) -> tuple[ShotParams, TableParams, tuple[float, float] | None]:
+ shot = ShotParams(
+ speed=req.speed,
+ angle=float(np.deg2rad(req.angle_deg)),
+ english_x=req.english_x,
+ english_y=req.english_y,
+ cue_elevation=float(np.deg2rad(req.cue_elevation_deg)),
+ )
+ table = TableParams(mu_slide=req.mu_slide, friction_noise_amp=req.friction_noise_amp)
+ obj = (req.obj_x, req.obj_y) if req.obj_x is not None and req.obj_y is not None else None
+ return shot, table, obj
+
+
@app.get("/health", response_model=HealthResponse)
def health() -> HealthResponse:
return HealthResponse(
status="ok",
- ml_loaded=_predictor.net is not None or _predictor.ort_session is not None,
- version="0.2.0",
+ ml_loaded=_predictor.ready,
+ backend=_predictor.backend,
+ version=API_VERSION,
)
@app.post("/predict")
def predict(req: ShotRequest) -> dict:
- import numpy as np
-
- shot = ShotParams(
- speed=req.speed,
- angle=float(np.deg2rad(req.angle_deg)),
- english_x=req.english_x,
- english_y=req.english_y,
- cue_elevation=float(np.deg2rad(req.cue_elevation_deg)),
- )
- table = TableParams(
- mu_slide=req.mu_slide,
- friction_noise_amp=req.friction_noise_amp,
- )
- obj = (req.obj_x, req.obj_y) if req.obj_x is not None and req.obj_y is not None else None
+ """Full numerical simulation, with the fast estimate alongside for comparison."""
+ shot, table, obj = _to_physics(req)
return _predictor.predict(
shot,
cue_pos=(req.cue_x, req.cue_y),
@@ -87,6 +105,18 @@ def predict(req: ShotRequest) -> dict:
)
+@app.post("/predict/fast")
+def predict_fast(req: ShotRequest) -> dict:
+ """Closed-form plus learned residual. No integration, sub-millisecond."""
+ shot, table, obj = _to_physics(req)
+ started = time.perf_counter()
+ result = _predictor.predict_fast(
+ shot, cue_pos=(req.cue_x, req.cue_y), obj_pos=obj, table=table, use_ml=req.use_ml
+ )
+ result["timing_ms"] = (time.perf_counter() - started) * 1000
+ return result
+
+
@app.get("/rack")
def rack(seed: int = 7) -> dict:
balls = make_full_rack(seed=seed)
diff --git a/src/cueai/ml/dataset.py b/src/cueai/ml/dataset.py
index e012613..630b85e 100644
--- a/src/cueai/ml/dataset.py
+++ b/src/cueai/ml/dataset.py
@@ -1,17 +1,34 @@
-"""Synthetic dataset generation from the physics simulator."""
+"""
+Synthetic dataset generation from the physics simulator.
+
+Each row pairs a shot with two things: the high-fidelity simulator outcome
+(the label) and the closed-form analytic prediction for the same shot (the
+baseline). The learned model predicts the gap between them, so its value is
+measured as error reduction over a real physical baseline rather than over
+nothing.
+
+Two properties are worth noting:
+
+* Table friction and restitution are resampled per shot (domain randomisation),
+ which is what makes the residual a function of the inputs rather than noise.
+* Each sample draws from its own seeded generator, so the dataset is identical
+ whether it is produced on 1 core or 32.
+"""
from __future__ import annotations
+import os
+from concurrent.futures import ProcessPoolExecutor
from pathlib import Path
import numpy as np
import pandas as pd
from tqdm import tqdm
+from cueai.physics.analytic import predict_endpoint
from cueai.physics.constants import ShotParams, TableParams
from cueai.physics.simulator import FEATURE_NAMES, Simulator, shot_feature_vector
-
TARGET_NAMES = [
"cue_end_x",
"cue_end_y",
@@ -21,6 +38,21 @@
"max_speed_cue",
]
+# Closed-form prediction for the same four endpoints, used as the residual base.
+BASELINE_NAMES = [
+ "baseline_cue_end_x",
+ "baseline_cue_end_y",
+ "baseline_obj_end_x",
+ "baseline_obj_end_y",
+]
+
+# Outcome complexity, recorded for analysis. These are results, not inputs, so
+# they must never be used as model features.
+DIAGNOSTIC_NAMES = ["n_cushion", "n_collision", "potted"]
+
+# Practical cue tip offset limit; beyond ~0.5R a real stroke miscues.
+MAX_TIP_OFFSET = 0.5
+
def _path_length(traj: np.ndarray) -> float:
if len(traj) < 2:
@@ -28,80 +60,104 @@ def _path_length(traj: np.ndarray) -> float:
return float(np.sum(np.linalg.norm(np.diff(traj, axis=0), axis=1)))
+def simulate_sample(index: int, seed: int = 42) -> dict[str, float]:
+ """
+ Simulate one random shot.
+
+ The generator is keyed on (seed, index) so the sample is reproducible in
+ isolation, independent of how the work is distributed across processes.
+ """
+ rng = np.random.default_rng([seed, index])
+
+ table = TableParams(
+ mu_slide=float(rng.uniform(0.15, 0.28)),
+ mu_roll=float(rng.uniform(0.006, 0.018)),
+ mu_spin=float(rng.uniform(0.03, 0.06)),
+ e_cushion=float(rng.uniform(0.75, 0.92)),
+ friction_noise_amp=float(rng.uniform(0.0, 0.04)),
+ )
+ sim = Simulator(table=table, dt=0.002, max_time=8.0)
+ length, width = table.length, table.width
+ radius = sim.ball_params.radius
+
+ cue_pos = np.array(
+ [rng.uniform(radius * 3, length * 0.4), rng.uniform(radius * 3, width - radius * 3)]
+ )
+ obj_pos = np.array(
+ [
+ rng.uniform(length * 0.45, length - radius * 3),
+ rng.uniform(radius * 3, width - radius * 3),
+ ]
+ )
+ if np.linalg.norm(cue_pos - obj_pos) < 4 * radius:
+ obj_pos[0] = min(length - 3 * radius, cue_pos[0] + 0.35)
+
+ shot = ShotParams(
+ speed=float(rng.uniform(0.5, 6.0)),
+ angle=float(rng.uniform(-np.pi, np.pi)),
+ english_x=float(rng.uniform(-MAX_TIP_OFFSET, MAX_TIP_OFFSET)),
+ english_y=float(rng.uniform(-MAX_TIP_OFFSET, MAX_TIP_OFFSET)),
+ cue_elevation=float(rng.uniform(0.0, 0.15)),
+ )
+ # Slightly over half the shots are aimed at the object ball, so the dataset
+ # is not dominated by shots that never make contact.
+ if rng.random() < 0.55:
+ delta = obj_pos - cue_pos
+ shot.angle = float(np.arctan2(delta[1], delta[0]) + rng.normal(0, 0.08))
+
+ result = sim.simulate_shot(
+ shot,
+ cue_pos=(float(cue_pos[0]), float(cue_pos[1])),
+ obj_pos=(float(obj_pos[0]), float(obj_pos[1])),
+ )
+
+ features = shot_feature_vector(shot, cue_pos, obj_pos, table)
+ cue_vel = result.velocities[0]
+ obj_end = result.endpoints.get(1, np.zeros(2))
+ targets = [
+ float(result.endpoints[0][0]),
+ float(result.endpoints[0][1]),
+ float(obj_end[0]),
+ float(obj_end[1]),
+ _path_length(result.trajectories[0]),
+ float(np.max(np.linalg.norm(cue_vel, axis=1))) if len(cue_vel) else 0.0,
+ ]
+
+ row = {name: float(value) for name, value in zip(FEATURE_NAMES, features)}
+ row.update(dict(zip(TARGET_NAMES, targets)))
+
+ # Closed-form baseline: analytic stopping point for the cue ball, and the
+ # object ball left where it started, since the analytic model has no notion
+ # of ball-ball contact.
+ baseline_cue = predict_endpoint(shot, cue_pos, table, radius=radius)
+ row["baseline_cue_end_x"] = float(baseline_cue[0])
+ row["baseline_cue_end_y"] = float(baseline_cue[1])
+ row["baseline_obj_end_x"] = float(obj_pos[0])
+ row["baseline_obj_end_y"] = float(obj_pos[1])
+
+ row["n_cushion"] = float(result.cushion_events)
+ row["n_collision"] = float(result.collision_events)
+ row["potted"] = float(sum(1 for was_potted in result.pocketed.values() if was_potted))
+ return row
+
+
def generate_dataset(
n_samples: int = 4000,
seed: int = 42,
out_csv: str | Path | None = "data/processed/shots.csv",
+ jobs: int | None = None,
) -> pd.DataFrame:
- rng = np.random.default_rng(seed)
- rows: list[dict] = []
-
- for _ in tqdm(range(n_samples), desc="Generating shots"):
- table = TableParams(
- mu_slide=float(rng.uniform(0.15, 0.28)),
- mu_roll=float(rng.uniform(0.006, 0.018)),
- mu_spin=float(rng.uniform(0.03, 0.06)),
- e_cushion=float(rng.uniform(0.75, 0.92)),
- friction_noise_amp=float(rng.uniform(0.0, 0.04)),
- )
- sim = Simulator(table=table, dt=0.002, max_time=8.0)
- L, W = table.length, table.width
- R = sim.ball_params.radius
- cue_pos = np.array(
- [rng.uniform(R * 3, L * 0.4), rng.uniform(R * 3, W - R * 3)],
- dtype=np.float64,
- )
- obj_pos = np.array(
- [rng.uniform(L * 0.45, L - R * 3), rng.uniform(R * 3, W - R * 3)],
- dtype=np.float64,
- )
- # Avoid overlapping start
- if np.linalg.norm(cue_pos - obj_pos) < 4 * R:
- obj_pos[0] = min(L - 3 * R, cue_pos[0] + 0.35)
-
- shot = ShotParams(
- speed=float(rng.uniform(0.5, 6.0)),
- angle=float(rng.uniform(-np.pi, np.pi)),
- english_x=float(rng.uniform(-0.6, 0.6)),
- english_y=float(rng.uniform(-0.6, 0.6)),
- cue_elevation=float(rng.uniform(0.0, 0.15)),
- )
-
- # Aim somewhat toward object ball sometimes
- if rng.random() < 0.55:
- delta = obj_pos - cue_pos
- shot.angle = float(np.arctan2(delta[1], delta[0]) + rng.normal(0, 0.08))
-
- result = sim.simulate_shot(
- shot,
- cue_pos=(float(cue_pos[0]), float(cue_pos[1])),
- obj_pos=(float(obj_pos[0]), float(obj_pos[1])),
- )
- x = shot_feature_vector(shot, cue_pos, obj_pos, table)
- cue_traj = result.trajectories[0]
- obj_traj = result.trajectories.get(1, np.zeros((1, 2)))
- cue_vel = result.velocities[0]
- y = np.array(
- [
- result.endpoints[0][0],
- result.endpoints[0][1],
- result.endpoints.get(1, np.zeros(2))[0],
- result.endpoints.get(1, np.zeros(2))[1],
- _path_length(cue_traj),
- float(np.max(np.linalg.norm(cue_vel, axis=1))) if len(cue_vel) else 0.0,
- ],
- dtype=np.float64,
- )
- # Physics baseline = endpoints without ML (same for residual training target)
- row = {n: float(v) for n, v in zip(FEATURE_NAMES, x)}
- for n, v in zip(TARGET_NAMES, y):
- row[n] = float(v)
- # Physics-only endpoint prediction used as baseline inside residual model
- row["phys_cue_end_x"] = float(result.endpoints[0][0])
- row["phys_cue_end_y"] = float(result.endpoints[0][1])
- row["phys_obj_end_x"] = float(result.endpoints.get(1, np.zeros(2))[0])
- row["phys_obj_end_y"] = float(result.endpoints.get(1, np.zeros(2))[1])
- rows.append(row)
+ """Simulate ``n_samples`` random shots, in parallel by default."""
+ jobs = jobs or min(os.cpu_count() or 1, 16)
+ indices = range(n_samples)
+ progress = {"total": n_samples, "desc": "Simulating shots", "unit": "shot"}
+
+ if jobs <= 1:
+ rows = [simulate_sample(i, seed) for i in tqdm(indices, **progress)]
+ else:
+ with ProcessPoolExecutor(max_workers=jobs) as pool:
+ futures = pool.map(simulate_sample, indices, [seed] * n_samples, chunksize=8)
+ rows = list(tqdm(futures, **progress))
df = pd.DataFrame(rows)
if out_csv is not None:
diff --git a/src/cueai/ml/features.py b/src/cueai/ml/features.py
new file mode 100644
index 0000000..bae88b2
--- /dev/null
+++ b/src/cueai/ml/features.py
@@ -0,0 +1,129 @@
+"""
+Feature construction, shared by training and serving.
+
+The raw shot parameters alone are a poor input for a residual model: to know
+where the closed-form prediction went wrong, the network would first have to
+rediscover cushion reflection geometry from a speed and an angle. So the
+features also carry what the closed-form solver already worked out — where it
+thinks the ball stops, how many cushions it expects, whether it expects a pot —
+plus the ghost-ball geometry that decides whether the object ball is touched at
+all.
+
+Two consequences worth stating:
+
+* The predicted cushion count tells the model when the outcome is chaotic, so it
+ can fall back to "trust the physics" instead of guessing.
+* Training and serving must build features identically. Both paths go through
+ :func:`build_features`, and ``tests/test_features.py`` pins them together.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+
+from cueai.physics.analytic import solve_free_ball
+from cueai.physics.constants import BallParams, ShotParams, TableParams
+from cueai.physics.simulator import FEATURE_NAMES as SHOT_FEATURE_NAMES
+
+BASELINE_FEATURE_NAMES = [
+ "base_cue_x", # closed-form resting position of the cue ball
+ "base_cue_y",
+ "base_cushions", # cushion contacts the closed-form solver expects
+ "base_potted", # whether it expects the cue ball to drop
+ "base_travel", # straight-line distance from cue start to that endpoint
+ "contact_along", # object ball distance projected onto the aim line
+ "contact_perp", # perpendicular miss distance of the aim line, signed
+ "will_contact", # aim line passes within one ball diameter, ahead of the cue
+]
+
+FEATURE_NAMES = SHOT_FEATURE_NAMES + BASELINE_FEATURE_NAMES
+
+
+def _spin_contact(shot: ShotParams, radius: float) -> np.ndarray:
+ """ω × (-Rẑ), the spin contribution to contact-point velocity."""
+ omega = shot.initial_omega(BallParams(radius=radius))
+ return np.array([-radius * omega[1], radius * omega[0]])
+
+
+def build_features(
+ shot: ShotParams,
+ cue_pos: np.ndarray | tuple[float, float],
+ obj_pos: np.ndarray | tuple[float, float] | None,
+ table: TableParams,
+ radius: float = 0.028575,
+) -> np.ndarray:
+ """Model input for a single shot. Order matches :data:`FEATURE_NAMES`."""
+ from cueai.physics.simulator import shot_feature_vector
+
+ cue = np.asarray(cue_pos, dtype=np.float64)
+ obj = np.asarray(obj_pos, dtype=np.float64) if obj_pos is not None else None
+
+ outcome = solve_free_ball(
+ cue,
+ shot.speed * np.array([np.cos(shot.angle), np.sin(shot.angle)]),
+ _spin_contact(shot, radius),
+ table,
+ radius,
+ )
+
+ aim = np.array([np.cos(shot.angle), np.sin(shot.angle)])
+ normal = np.array([-aim[1], aim[0]])
+ if obj is None:
+ along, perp = 0.0, 0.0
+ else:
+ offset = obj - cue
+ along, perp = float(offset @ aim), float(offset @ normal)
+ will_contact = float(along > 0 and abs(perp) < 2 * radius)
+
+ return np.concatenate(
+ [
+ shot_feature_vector(shot, cue, obj, table),
+ [
+ outcome.position[0],
+ outcome.position[1],
+ float(outcome.cushions),
+ float(outcome.potted),
+ float(np.linalg.norm(outcome.position - cue)),
+ along,
+ perp,
+ will_contact,
+ ],
+ ]
+ )
+
+
+def build_feature_frame(df: pd.DataFrame, radius: float = 0.028575) -> np.ndarray:
+ """
+ Rebuild the feature matrix from a generated dataset.
+
+ The stored CSV keeps raw shot and table parameters rather than derived
+ features, so feature engineering can change without re-running 20,000
+ simulations.
+ """
+ rows = []
+ for record in df.to_dict("records"):
+ shot = ShotParams(
+ speed=float(record["speed"]),
+ angle=float(record["angle"]),
+ english_x=float(record["english_x"]),
+ english_y=float(record["english_y"]),
+ cue_elevation=float(record["cue_elevation"]),
+ )
+ table = TableParams(
+ mu_slide=float(record["mu_slide"]),
+ mu_roll=float(record["mu_roll"]),
+ mu_spin=float(record["mu_spin"]),
+ e_cushion=float(record["e_cushion"]),
+ friction_noise_amp=float(record["friction_noise_amp"]),
+ )
+ rows.append(
+ build_features(
+ shot,
+ (float(record["cue_x"]), float(record["cue_y"])),
+ (float(record["obj_x"]), float(record["obj_y"])),
+ table,
+ radius,
+ )
+ )
+ return np.asarray(rows, dtype=np.float64)
diff --git a/src/cueai/ml/infer.py b/src/cueai/ml/infer.py
index 8bcb018..59e000b 100644
--- a/src/cueai/ml/infer.py
+++ b/src/cueai/ml/infer.py
@@ -1,19 +1,36 @@
-"""Inference helpers: physics + ONNX / Torch residual fusion."""
+"""
+Inference: the closed-form baseline, the learned residual, and the simulator.
+
+Two prediction paths are exposed deliberately, because they trade accuracy for
+latency by four orders of magnitude:
+
+``predict_fast``
+ Closed-form solution plus the CueNet residual. Sub-millisecond, endpoints
+ only, suitable for search or for a real-time aiming aid.
+``predict``
+ Full numerical simulation for the whole rack, which is what the desktop UI
+ animates, alongside the fast prediction and the gap between them.
+"""
from __future__ import annotations
+import time
from pathlib import Path
import numpy as np
import torch
+from cueai.ml.features import FEATURE_NAMES, build_features
from cueai.ml.model import CueNet
+from cueai.physics.analytic import predict_endpoint
from cueai.physics.ball import Ball
from cueai.physics.constants import ShotParams, TableParams
-from cueai.physics.simulator import FEATURE_NAMES, Simulator, shot_feature_vector
+from cueai.physics.simulator import Simulator
class TrajectoryPredictor:
+ """Loads whichever CueNet artefacts are present and serves predictions."""
+
def __init__(self, model_dir: str | Path = "models"):
self.model_dir = Path(model_dir)
self.sim = Simulator(dt=0.001, max_time=14.0, collision_passes=20)
@@ -23,41 +40,106 @@ def __init__(self, model_dir: str | Path = "models"):
self.ort_session = None
self._load()
+ # ------------------------------------------------------------------ loading
+
def _load(self) -> None:
- ckpt = self.model_dir / "cuenet.pt"
- onnx = self.model_dir / "cuenet.onnx"
- if ckpt.exists():
- data = torch.load(ckpt, map_location="cpu", weights_only=False)
- self.net = CueNet(in_dim=int(data["in_dim"]))
+ checkpoint = self.model_dir / "cuenet.pt"
+ onnx_path = self.model_dir / "cuenet.onnx"
+ if checkpoint.exists():
+ data = torch.load(checkpoint, map_location="cpu", weights_only=False)
+ self.net = CueNet(in_dim=int(data["in_dim"]), hidden=int(data.get("hidden", 256)))
self.net.load_state_dict(data["state_dict"])
self.net.eval()
self.mean = np.asarray(data["scaler_mean"], dtype=np.float32)
self.scale = np.asarray(data["scaler_scale"], dtype=np.float32)
- if onnx.exists():
+ if onnx_path.exists():
try:
import onnxruntime as ort
self.ort_session = ort.InferenceSession(
- str(onnx), providers=["CPUExecutionProvider"]
+ str(onnx_path), providers=["CPUExecutionProvider"]
)
- except Exception:
+ except Exception: # pragma: no cover - optional runtime
self.ort_session = None
- def _scale(self, x: np.ndarray) -> np.ndarray:
+ @property
+ def ready(self) -> bool:
+ """True when a trained residual model is available."""
+ return self.net is not None or self.ort_session is not None
+
+ @property
+ def backend(self) -> str:
+ if self.ort_session is not None:
+ return "onnx"
+ return "torch" if self.net is not None else "none"
+
+ # --------------------------------------------------------------- prediction
+
+ def feature_vector(
+ self,
+ shot: ShotParams,
+ cue_pos: tuple[float, float] | np.ndarray,
+ obj_pos: tuple[float, float] | np.ndarray | None,
+ table: TableParams,
+ ) -> np.ndarray:
+ """Identical construction to training; see :mod:`cueai.ml.features`."""
+ return build_features(
+ shot, cue_pos, obj_pos, table, radius=self.sim.ball_params.radius
+ )
+
+ def _standardise(self, features: np.ndarray) -> np.ndarray:
if self.mean is None or self.scale is None:
- return x.astype(np.float32)
- return ((x - self.mean) / np.clip(self.scale, 1e-8, None)).astype(np.float32)
+ return features.astype(np.float32)
+ return ((features - self.mean) / np.clip(self.scale, 1e-8, None)).astype(np.float32)
- def _residual(self, x: np.ndarray) -> np.ndarray:
- xs = self._scale(x)[None, :]
+ def residual_batch(self, features: np.ndarray) -> np.ndarray:
+ """Endpoint corrections for a batch of raw feature rows, shape (N, 4)."""
+ standardised = self._standardise(np.atleast_2d(features))
if self.ort_session is not None:
- out = self.ort_session.run(None, {"features": xs})[0]
- return np.asarray(out[0], dtype=np.float64)
+ return np.asarray(
+ self.ort_session.run(None, {"features": standardised})[0], dtype=np.float64
+ )
if self.net is not None:
with torch.no_grad():
- out = self.net(torch.from_numpy(xs)).numpy()[0]
- return out.astype(np.float64)
- return np.zeros(4, dtype=np.float64)
+ return self.net(torch.from_numpy(standardised)).numpy().astype(np.float64)
+ return np.zeros((len(standardised), 4), dtype=np.float64)
+
+ def predict_fast(
+ self,
+ shot: ShotParams,
+ cue_pos: tuple[float, float],
+ obj_pos: tuple[float, float] | None = None,
+ table: TableParams | None = None,
+ use_ml: bool = True,
+ ) -> dict:
+ """Closed-form endpoints plus the learned residual. No integration."""
+ table = table or self.sim.table
+ baseline_cue = predict_endpoint(shot, cue_pos, table, radius=self.sim.ball_params.radius)
+ baseline_obj = np.asarray(obj_pos, dtype=np.float64) if obj_pos is not None else None
+ baseline = np.concatenate(
+ [baseline_cue, baseline_obj if baseline_obj is not None else np.zeros(2)]
+ )
+
+ features = self.feature_vector(shot, cue_pos, obj_pos, table)
+ residual = (
+ self.residual_batch(features)[0] if (use_ml and self.ready) else np.zeros(4)
+ )
+ corrected = baseline + residual
+ # With no object ball there is nothing for the object channel to refer to,
+ # so it is reported as absent rather than as a position near the origin.
+ has_object = baseline_obj is not None
+ return {
+ "baseline_endpoints": {
+ "cue": baseline[:2].tolist(),
+ "object": baseline[2:].tolist() if has_object else None,
+ },
+ "residual": residual.tolist(),
+ "endpoints": {
+ "cue": corrected[:2].tolist(),
+ "object": corrected[2:].tolist() if has_object else None,
+ },
+ "backend": self.backend if use_ml else "closed_form",
+ }
def predict(
self,
@@ -70,8 +152,16 @@ def predict(
seed: int | None = 7,
balls: list[Ball] | None = None,
) -> dict:
+ """
+ Full simulation plus the fast prediction, so the two can be compared.
+
+ The returned ``agreement`` block is what the UI and the API surface: how
+ far the sub-millisecond estimate lands from the simulated outcome.
+ """
if table is not None:
self.sim.table = table
+
+ sim_start = time.perf_counter()
result = self.sim.simulate_shot(
shot,
cue_pos=cue_pos,
@@ -80,39 +170,42 @@ def predict(
seed=seed,
balls=balls,
)
- cue_p = np.array(cue_pos, dtype=np.float64)
- # Use 8-ball as reference object for ML residual (legacy head)
- eight = result.endpoints.get(8, result.endpoints.get(1, np.zeros(2)))
- obj_p = np.asarray(eight, dtype=np.float64)
- feats = shot_feature_vector(shot, cue_p, obj_p, self.sim.table)
- phys = np.array(
- [
- result.endpoints.get(0, np.zeros(2))[0],
- result.endpoints.get(0, np.zeros(2))[1],
- float(obj_p[0]),
- float(obj_p[1]),
- ]
+ sim_ms = (time.perf_counter() - sim_start) * 1000
+
+ # With a full rack the 8-ball stands in for "the object ball" so that the
+ # single-object-ball model trained on two-ball shots still has a referent.
+ reference = result.endpoints.get(8, result.endpoints.get(1, np.zeros(2)))
+ reference_pos = (float(reference[0]), float(reference[1]))
+
+ fast_start = time.perf_counter()
+ fast = self.predict_fast(
+ shot, cue_pos, obj_pos=obj_pos or reference_pos, table=self.sim.table, use_ml=use_ml
)
- residual = self._residual(feats) if use_ml else np.zeros(4)
- corrected = phys + residual
+ fast_ms = (time.perf_counter() - fast_start) * 1000
+
+ simulated_cue = np.asarray(result.endpoints.get(0, np.zeros(2)), dtype=np.float64)
+ cue_gap = float(np.linalg.norm(simulated_cue - np.asarray(fast["endpoints"]["cue"])))
+
return {
- "features": {k: float(v) for k, v in zip(FEATURE_NAMES, feats)},
- "physics_endpoints": {
- "cue": phys[:2].tolist(),
- "object": phys[2:].tolist(),
- },
- "ml_residual": residual.tolist(),
- "fused_endpoints": {
- "cue": corrected[:2].tolist(),
- "object": corrected[2:].tolist(),
+ "features": dict(
+ zip(
+ FEATURE_NAMES,
+ self.feature_vector(shot, cue_pos, reference_pos, self.sim.table),
+ )
+ ),
+ "simulated_endpoints": {
+ "cue": simulated_cue.tolist(),
+ "object": list(reference_pos),
},
+ "fast_prediction": fast,
+ "agreement": {"cue_gap_m": cue_gap},
+ "timing_ms": {"simulator": sim_ms, "fast": fast_ms},
"endpoints": {str(k): v.tolist() for k, v in result.endpoints.items()},
- "trajectory": {
- str(k): v.tolist() for k, v in result.trajectories.items()
- },
+ "trajectory": {str(k): v.tolist() for k, v in result.trajectories.items()},
"ball_meta": {str(k): v for k, v in result.ball_meta.items()},
"pocketed": {str(k): bool(v) for k, v in result.pocketed.items()},
"collisions": result.collision_events,
+ "cushions": result.cushion_events,
"times": result.times.tolist(),
- "ml_loaded": self.net is not None or self.ort_session is not None,
+ "ml_loaded": self.ready,
}
diff --git a/src/cueai/ml/model.py b/src/cueai/ml/model.py
index a8b7f55..c4e2c64 100644
--- a/src/cueai/ml/model.py
+++ b/src/cueai/ml/model.py
@@ -16,6 +16,11 @@ class CueNet(nn.Module):
def __init__(self, in_dim: int = 18, hidden: int = 128, out_dim: int = 4):
super().__init__()
+ # Zero-initialised head: training starts from "trust the physics exactly"
+ # and only moves away where the data says the physics is incomplete.
+ head = nn.Linear(hidden // 2, out_dim)
+ nn.init.zeros_(head.weight)
+ nn.init.zeros_(head.bias)
self.net = nn.Sequential(
nn.Linear(in_dim, hidden),
nn.SiLU(),
@@ -25,11 +30,8 @@ def __init__(self, in_dim: int = 18, hidden: int = 128, out_dim: int = 4):
nn.Dropout(0.1),
nn.Linear(hidden, hidden // 2),
nn.SiLU(),
- nn.Linear(hidden // 2, out_dim),
+ head,
)
- # Start near identity residual (physics-only)
- nn.init.zeros_(self.net[-1].weight)
- nn.init.zeros_(self.net[-1].bias)
def forward(self, x: torch.Tensor) -> torch.Tensor:
return self.net(x)
diff --git a/src/cueai/ml/train.py b/src/cueai/ml/train.py
index 38f8745..7548277 100644
--- a/src/cueai/ml/train.py
+++ b/src/cueai/ml/train.py
@@ -1,9 +1,22 @@
-"""Train CueNet residual model + export ONNX + sklearn baseline."""
+"""
+Train the CueNet residual model, export ONNX, and score it against baselines.
+
+The learned quantity is the gap between a closed-form prediction and the
+high-fidelity simulator, so every reported number is an error reduction over a
+physical baseline. Three models are compared on the same held-out split:
+
+ analytic closed-form, no fitting (the baseline to beat)
+ gbm gradient boosting on raw features (does ML help at all?)
+ cuenet MLP predicting the analytic residual
+
+Results land in ``models/metrics.json``, which is what the README quotes.
+"""
from __future__ import annotations
import argparse
import json
+import platform
from pathlib import Path
import joblib
@@ -11,150 +24,284 @@
import pandas as pd
import torch
import torch.nn as nn
-from sklearn.ensemble import GradientBoostingRegressor
-from sklearn.multioutput import MultiOutputRegressor
+from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.model_selection import train_test_split
+from sklearn.multioutput import MultiOutputRegressor
from sklearn.preprocessing import StandardScaler
from torch.utils.data import DataLoader, TensorDataset
-from cueai.ml.dataset import FEATURE_NAMES, TARGET_NAMES, generate_dataset
+from cueai.ml.dataset import BASELINE_NAMES, TARGET_NAMES, generate_dataset
+from cueai.ml.features import BASELINE_FEATURE_NAMES, FEATURE_NAMES, build_feature_frame
from cueai.ml.model import CueNet
+# Endpoint targets, in metres: (cue_x, cue_y, obj_x, obj_y)
+ENDPOINT_TARGETS = TARGET_NAMES[:4]
+TEST_SIZE = 0.2
+SPLIT_SEED = 0
+
+
+def endpoint_errors(pred: np.ndarray, truth: np.ndarray) -> dict[str, float]:
+ """
+ Error summary in millimetres.
+
+ ``euclidean`` is the mean distance between predicted and true resting
+ positions, averaged over the cue and object ball, which is the number a
+ player or a downstream planner actually cares about.
+ """
+ per_axis = np.abs(pred - truth)
+ cue_dist = np.linalg.norm(pred[:, :2] - truth[:, :2], axis=1)
+ obj_dist = np.linalg.norm(pred[:, 2:] - truth[:, 2:], axis=1)
+ both = np.concatenate([cue_dist, obj_dist])
+ ss_res = float(np.sum((pred - truth) ** 2))
+ ss_tot = float(np.sum((truth - truth.mean(axis=0)) ** 2))
+ return {
+ "mae_mm": float(per_axis.mean() * 1000),
+ "euclidean_mm": float(both.mean() * 1000),
+ "p95_mm": float(np.percentile(both, 95) * 1000),
+ "cue_mm": float(cue_dist.mean() * 1000),
+ "obj_mm": float(obj_dist.mean() * 1000),
+ "r2": float(1.0 - ss_res / ss_tot) if ss_tot > 0 else float("nan"),
+ }
+
+
+def split_indices(n_rows: int) -> tuple[np.ndarray, np.ndarray]:
+ """One split shared by every model, so the comparison is apples to apples."""
+ return train_test_split(
+ np.arange(n_rows), test_size=TEST_SIZE, random_state=SPLIT_SEED
+ )
-PHYS_COLS = ["phys_cue_end_x", "phys_cue_end_y", "phys_obj_end_x", "phys_obj_end_y"]
-RESIDUAL_TARGETS = TARGET_NAMES[:4] # cue/obj endpoints
+def stratify_by_contacts(
+ predictions: dict[str, np.ndarray], truth: np.ndarray, n_cushion: np.ndarray
+) -> list[dict]:
+ """
+ Break accuracy down by how many cushions the cue ball touched.
-def train_sklearn_baseline(df: pd.DataFrame, model_dir: Path) -> dict:
- X = df[FEATURE_NAMES].values
- y = df[RESIDUAL_TARGETS].values
- X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=0)
+ Resting position is a smooth function of the shot until the ball starts
+ ricocheting; past a couple of cushion contacts the outcome is chaotic and no
+ model should be expected to predict it. Reporting the breakdown states the
+ limit of the approach instead of hiding it in an average.
+ """
+ buckets = [("0", 0, 0), ("1", 1, 1), ("2", 2, 2), ("3+", 3, 10_000)]
+ rows = []
+ for label, low, high in buckets:
+ mask = (n_cushion >= low) & (n_cushion <= high)
+ if not mask.any():
+ continue
+ row: dict = {"cushion_contacts": label, "n": int(mask.sum())}
+ for name, pred in predictions.items():
+ row[name] = endpoint_errors(pred[mask], truth[mask])["euclidean_mm"]
+ # How much of the true spread the model is willing to reproduce. A ratio
+ # well below 1 means it is hedging toward the middle of the table, which
+ # is the right move when the outcome is genuinely unpredictable.
+ row["cuenet_spread_ratio"] = float(
+ predictions["cuenet"][mask].std() / max(truth[mask].std(), 1e-9)
+ )
+ rows.append(row)
+ return rows
+
+
+def stratify_by_object_contact(
+ predictions: dict[str, np.ndarray], truth: np.ndarray, n_collision: np.ndarray
+) -> list[dict]:
+ """
+ Split accuracy by whether the cue ball ever touched the object ball.
+
+ This matters for reading the headline number honestly. Sampling aims roughly
+ half the shots at the object ball and misses often, so most shots leave it
+ where it started — and since the closed-form baseline predicts exactly that,
+ the object-ball half of the reported error is trivially satisfied for them.
+ Averaging the two subsets together would flatter every model equally; the
+ split shows what each one does when a collision actually has to be modelled.
+ """
+ rows = []
+ for label, mask in (("no", n_collision == 0), ("yes", n_collision > 0)):
+ if not mask.any():
+ continue
+ row: dict = {
+ "object_ball_contact": label,
+ "n": int(mask.sum()),
+ "share_pct": round(100.0 * float(mask.mean()), 1),
+ }
+ for name, pred in predictions.items():
+ errors = endpoint_errors(pred[mask], truth[mask])
+ row[f"{name}_cue_mm"] = errors["cue_mm"]
+ row[f"{name}_obj_mm"] = errors["obj_mm"]
+ rows.append(row)
+ return rows
+
+
+def risk_coverage(
+ predictions: dict[str, np.ndarray], truth: np.ndarray, expected_cushions: np.ndarray
+) -> list[dict]:
+ """
+ Accuracy as a function of how much of the shot space the model answers for.
+
+ :func:`stratify_by_contacts` slices by what the simulator *did*, which is only
+ knowable after paying for the simulation. This slices by what the closed-form
+ solver *expects* before anything is run, so it is usable as a gate: answer
+ only when the solver predicts at most ``k`` cushion contacts, and defer the
+ rest to the simulator. The gate is free, because the expected cushion count
+ is already computed as part of producing the prediction.
+ """
+ rows = []
+ for threshold in (0, 1, 2, 3, None):
+ mask = (
+ np.ones(len(truth), dtype=bool)
+ if threshold is None
+ else expected_cushions <= threshold
+ )
+ if not mask.any():
+ continue
+ row: dict = {
+ "expected_cushions_at_most": "all" if threshold is None else str(threshold),
+ "n": int(mask.sum()),
+ "coverage_pct": round(100.0 * float(mask.mean()), 1),
+ }
+ for name, pred in predictions.items():
+ row[name] = endpoint_errors(pred[mask], truth[mask])["euclidean_mm"]
+ rows.append(row)
+ return rows
+
+
+def train_gbm(
+ x_train: np.ndarray,
+ y_train: np.ndarray,
+ x_test: np.ndarray,
+ y_test: np.ndarray,
+ model_dir: Path,
+) -> tuple[np.ndarray, dict[str, float]]:
scaler = StandardScaler()
- X_tr_s = scaler.fit_transform(X_tr)
- X_te_s = scaler.transform(X_te)
model = MultiOutputRegressor(
- GradientBoostingRegressor(n_estimators=120, max_depth=4, random_state=0)
+ HistGradientBoostingRegressor(max_iter=300, learning_rate=0.1, random_state=0)
)
- model.fit(X_tr_s, y_tr)
- pred = model.predict(X_te_s)
- mae = float(np.mean(np.abs(pred - y_te)))
- joblib.dump({"model": model, "scaler": scaler}, model_dir / "sklearn_baseline.joblib")
- return {"sklearn_mae_m": mae}
+ model.fit(scaler.fit_transform(x_train), y_train)
+ pred = model.predict(scaler.transform(x_test))
+ joblib.dump({"model": model, "scaler": scaler}, model_dir / "gbm_baseline.joblib")
+ return pred, endpoint_errors(pred, y_test)
-def train_torch(
- df: pd.DataFrame,
+def train_cuenet(
+ x_train: np.ndarray,
+ base_train: np.ndarray,
+ y_train: np.ndarray,
+ x_test: np.ndarray,
+ base_test: np.ndarray,
+ y_test: np.ndarray,
model_dir: Path,
- epochs: int = 40,
- batch_size: int = 128,
- lr: float = 1e-3,
-) -> dict:
- device = torch.device("cpu")
- X = df[FEATURE_NAMES].values.astype(np.float32)
- phys = df[PHYS_COLS].values.astype(np.float32)
- y = df[RESIDUAL_TARGETS].values.astype(np.float32)
- # Train residual = y - phys (ideally ~0; noise/table variance creates learnable signal
- # when we inject measurement noise)
- noise = np.random.default_rng(1).normal(0, 0.008, size=y.shape).astype(np.float32)
- y_obs = y + noise
- residual = y_obs - phys
-
- X_tr, X_te, p_tr, p_te, r_tr, r_te, y_tr, y_te = train_test_split(
- X, phys, residual, y_obs, test_size=0.2, random_state=0
+ epochs: int,
+ batch_size: int,
+ lr: float,
+ seed: int,
+ hidden: int = 256,
+ save: bool = True,
+) -> tuple[np.ndarray, dict[str, float], list[float]]:
+ """
+ Fit the residual model, selecting the epoch on a validation split.
+
+ The test split is scored exactly once, after training, so no model choice is
+ informed by it. ``save=False`` fits a throwaway model for the feature
+ ablation, so the diagnostic cannot overwrite the served artefacts.
+ """
+ torch.manual_seed(seed)
+ fit_x, val_x, fit_base, val_base, fit_y, val_y = train_test_split(
+ x_train, base_train, y_train, test_size=0.1, random_state=SPLIT_SEED
)
+
scaler = StandardScaler()
- X_tr_s = scaler.fit_transform(X_tr).astype(np.float32)
- X_te_s = scaler.transform(X_te).astype(np.float32)
+ fit_x_s = scaler.fit_transform(fit_x).astype(np.float32)
+ val_x_s = scaler.transform(val_x).astype(np.float32)
+ x_test_s = scaler.transform(x_test).astype(np.float32)
- ds = TensorDataset(
- torch.from_numpy(X_tr_s),
- torch.from_numpy(p_tr),
- torch.from_numpy(r_tr),
+ # The network only has to explain what the closed-form model misses: cushion
+ # friction and spin transfer, ball-ball contact, throw, and cloth variation.
+ loader = DataLoader(
+ TensorDataset(torch.from_numpy(fit_x_s), torch.from_numpy(fit_y - fit_base)),
+ batch_size=batch_size,
+ shuffle=True,
)
- loader = DataLoader(ds, batch_size=batch_size, shuffle=True)
-
- net = CueNet(in_dim=X.shape[1]).to(device)
+ net = CueNet(in_dim=x_train.shape[1], hidden=hidden)
opt = torch.optim.AdamW(net.parameters(), lr=lr, weight_decay=1e-4)
- loss_fn = nn.SmoothL1Loss()
-
- history = []
- net.train()
- for epoch in range(epochs):
- total = 0.0
- n = 0
- for xb, _pb, rb in loader:
- xb, rb = xb.to(device), rb.to(device)
- pred_r = net(xb)
- loss = loss_fn(pred_r, rb)
+ sched = torch.optim.lr_scheduler.CosineAnnealingLR(opt, T_max=max(epochs, 1))
+ loss_fn = nn.SmoothL1Loss(beta=0.05)
+
+ history: list[float] = []
+ best_error, best_state = np.inf, None
+ val_tensor = torch.from_numpy(val_x_s)
+ for _ in range(epochs):
+ net.train()
+ running, count = 0.0, 0
+ for features, residual in loader:
+ loss = loss_fn(net(features), residual)
opt.zero_grad()
loss.backward()
opt.step()
- total += float(loss.item()) * len(xb)
- n += len(xb)
- history.append(total / max(n, 1))
+ running += float(loss.item()) * len(features)
+ count += len(features)
+ sched.step()
+ history.append(running / max(count, 1))
+
+ net.eval()
+ with torch.no_grad():
+ val_pred = val_base + net(val_tensor).numpy()
+ val_error = endpoint_errors(val_pred, val_y)["euclidean_mm"]
+ if val_error < best_error:
+ best_error = val_error
+ best_state = {k: v.clone() for k, v in net.state_dict().items()}
+ if best_state is not None:
+ net.load_state_dict(best_state)
net.eval()
with torch.no_grad():
- pred_r = net(torch.from_numpy(X_te_s)).numpy()
- pred_y = p_te + pred_r
- mae = float(np.mean(np.abs(pred_y - y_te)))
- phys_mae = float(np.mean(np.abs(p_te - y_te)))
-
- ckpt = model_dir / "cuenet.pt"
- torch.save(
- {
- "state_dict": net.state_dict(),
- "scaler_mean": scaler.mean_,
- "scaler_scale": scaler.scale_,
- "feature_names": FEATURE_NAMES,
- "in_dim": X.shape[1],
- },
- ckpt,
- )
+ pred = base_test + net(torch.from_numpy(x_test_s)).numpy()
+
+ if save:
+ torch.save(
+ {
+ "state_dict": net.state_dict(),
+ "scaler_mean": scaler.mean_,
+ "scaler_scale": scaler.scale_,
+ "feature_names": FEATURE_NAMES,
+ "in_dim": x_train.shape[1],
+ "hidden": hidden,
+ },
+ model_dir / "cuenet.pt",
+ )
+ _export_onnx(net, x_train.shape[1], model_dir / "cuenet.onnx")
+ return pred, endpoint_errors(pred, y_test), history
- # ONNX export (best-effort; Torch 2.x may need onnxscript)
- onnx_path = model_dir / "cuenet.onnx"
- dummy = torch.randn(1, X.shape[1], dtype=torch.float32)
+
+def _export_onnx(net: CueNet, in_dim: int, path: Path) -> bool:
+ """Export for runtime-agnostic serving. Non-fatal: the .pt checkpoint stands alone."""
try:
torch.onnx.export(
net.cpu(),
- dummy,
- str(onnx_path),
+ (torch.randn(1, in_dim, dtype=torch.float32),),
+ str(path),
input_names=["features"],
output_names=["residual"],
dynamic_axes={"features": {0: "batch"}, "residual": {0: "batch"}},
opset_version=17,
dynamo=False,
)
- except Exception as exc: # noqa: BLE001
+ return True
+ except Exception as exc:
print(f"ONNX export skipped: {exc}")
- onnx_path = model_dir / "cuenet.onnx"
- if not onnx_path.exists():
- onnx_path = Path("")
+ return False
- meta = {
- "torch_mae_m": mae,
- "physics_only_mae_m": phys_mae,
- "improvement_pct": float(100 * (phys_mae - mae) / max(phys_mae, 1e-9)),
- "epochs": epochs,
- "n_samples": len(df),
- "history_last": history[-1] if history else None,
- "onnx": str(onnx_path),
- "checkpoint": str(ckpt),
- }
- (model_dir / "metrics.json").write_text(json.dumps(meta, indent=2))
- return meta
-
-
-def main(argv: list[str] | None = None) -> None:
- p = argparse.ArgumentParser(description="Train CueAI models")
- p.add_argument("--n-samples", type=int, default=3000)
- p.add_argument("--epochs", type=int, default=40)
- p.add_argument("--data", type=str, default="data/processed/shots.csv")
- p.add_argument("--model-dir", type=str, default="models")
- p.add_argument("--skip-generate", action="store_true")
- args = p.parse_args(argv)
+def main(argv: list[str] | None = None) -> dict:
+ parser = argparse.ArgumentParser(description="Train CueAI models")
+ parser.add_argument("--n-samples", type=int, default=4000)
+ parser.add_argument("--epochs", type=int, default=120)
+ parser.add_argument("--batch-size", type=int, default=128)
+ parser.add_argument("--hidden", type=int, default=256)
+ parser.add_argument("--lr", type=float, default=2e-3)
+ parser.add_argument("--seed", type=int, default=0)
+ parser.add_argument("--data", type=str, default="data/processed/shots.csv")
+ parser.add_argument("--model-dir", type=str, default="models")
+ parser.add_argument("--skip-generate", action="store_true")
+ args = parser.parse_args(argv)
model_dir = Path(args.model_dir)
model_dir.mkdir(parents=True, exist_ok=True)
@@ -165,10 +312,102 @@ def main(argv: list[str] | None = None) -> None:
else:
df = generate_dataset(n_samples=args.n_samples, out_csv=data_path)
- sk = train_sklearn_baseline(df, model_dir)
- torch_meta = train_torch(df, model_dir, epochs=args.epochs)
- summary = {**sk, **torch_meta}
- print(json.dumps(summary, indent=2))
+ features = build_feature_frame(df).astype(np.float32)
+ baselines = df[BASELINE_NAMES].to_numpy(np.float32)
+ targets = df[ENDPOINT_TARGETS].to_numpy(np.float32)
+ train_idx, test_idx = split_indices(len(df))
+ x_train, x_test = features[train_idx], features[test_idx]
+ base_train, base_test = baselines[train_idx], baselines[test_idx]
+ y_train, y_test = targets[train_idx], targets[test_idx]
+
+ analytic = endpoint_errors(base_test, y_test)
+ gbm_pred, gbm = train_gbm(x_train, y_train, x_test, y_test, model_dir)
+ cuenet_pred, cuenet, history = train_cuenet(
+ x_train,
+ base_train,
+ y_train,
+ x_test,
+ base_test,
+ y_test,
+ model_dir,
+ epochs=args.epochs,
+ batch_size=args.batch_size,
+ lr=args.lr,
+ seed=args.seed,
+ hidden=args.hidden,
+ )
+
+ # Ablation: the same network, epochs and split, but without the closed-form
+ # solver's conclusions among its inputs. This is the measurement behind the
+ # claim that the features, not the architecture, are what let the residual
+ # model beat the physics on the shots physics nearly solves.
+ n_shot_features = len(FEATURE_NAMES) - len(BASELINE_FEATURE_NAMES)
+ ablation_pred, ablation, _ = train_cuenet(
+ x_train[:, :n_shot_features],
+ base_train,
+ y_train,
+ x_test[:, :n_shot_features],
+ base_test,
+ y_test,
+ model_dir,
+ epochs=args.epochs,
+ batch_size=args.batch_size,
+ lr=args.lr,
+ seed=args.seed,
+ hidden=args.hidden,
+ save=False,
+ )
+
+ predictions = {"analytic": base_test, "gbm": gbm_pred, "cuenet": cuenet_pred}
+ direct = df["n_cushion"].to_numpy()[test_idx] == 0
+ metrics = {
+ "n_samples": len(df),
+ "n_test": len(y_test),
+ "epochs": args.epochs,
+ "hidden": args.hidden,
+ "final_train_loss": history[-1] if history else None,
+ "models": {"analytic": analytic, "gbm": gbm, "cuenet": cuenet},
+ "by_cushion_contacts": stratify_by_contacts(
+ predictions, y_test, df["n_cushion"].to_numpy()[test_idx]
+ ),
+ "by_object_ball_contact": stratify_by_object_contact(
+ predictions, y_test, df["n_collision"].to_numpy()[test_idx]
+ ),
+ "risk_coverage": risk_coverage(
+ predictions, y_test, x_test[:, FEATURE_NAMES.index("base_cushions")]
+ ),
+ "feature_ablation": {
+ "note": (
+ "identical architecture, epochs, seed and split; the ablated model "
+ "sees only raw shot parameters, not the closed-form solver's output"
+ ),
+ "shot_features_only_mm": ablation["euclidean_mm"],
+ "with_closed_form_features_mm": cuenet["euclidean_mm"],
+ "shot_features_only_direct_mm": endpoint_errors(
+ ablation_pred[direct], y_test[direct]
+ )["euclidean_mm"],
+ "with_closed_form_features_direct_mm": endpoint_errors(
+ cuenet_pred[direct], y_test[direct]
+ )["euclidean_mm"],
+ "analytic_direct_mm": endpoint_errors(base_test[direct], y_test[direct])[
+ "euclidean_mm"
+ ],
+ },
+ "error_reduction_vs_analytic_pct": round(
+ 100 * (1 - cuenet["euclidean_mm"] / analytic["euclidean_mm"]), 1
+ ),
+ "error_reduction_vs_gbm_pct": round(
+ 100 * (1 - cuenet["euclidean_mm"] / gbm["euclidean_mm"]), 1
+ ),
+ "environment": {
+ "python": platform.python_version(),
+ "torch": torch.__version__,
+ "platform": platform.platform(),
+ },
+ }
+ (model_dir / "metrics.json").write_text(json.dumps(metrics, indent=2) + "\n")
+ print(json.dumps(metrics, indent=2))
+ return metrics
if __name__ == "__main__":
diff --git a/src/cueai/physics/__init__.py b/src/cueai/physics/__init__.py
index 3021d9c..2807080 100644
--- a/src/cueai/physics/__init__.py
+++ b/src/cueai/physics/__init__.py
@@ -1,15 +1,18 @@
"""Physics package."""
+from cueai.physics.analytic import predict_endpoint, straight_shot
from cueai.physics.constants import BallParams, ShotParams, TableParams
from cueai.physics.rack import make_full_rack
-from cueai.physics.simulator import Simulator, SimResult, shot_feature_vector
+from cueai.physics.simulator import SimResult, Simulator, shot_feature_vector
__all__ = [
"BallParams",
"ShotParams",
- "TableParams",
- "Simulator",
"SimResult",
- "shot_feature_vector",
+ "Simulator",
+ "TableParams",
"make_full_rack",
+ "predict_endpoint",
+ "shot_feature_vector",
+ "straight_shot",
]
diff --git a/src/cueai/physics/analytic.py b/src/cueai/physics/analytic.py
new file mode 100644
index 0000000..c21b99d
--- /dev/null
+++ b/src/cueai/physics/analytic.py
@@ -0,0 +1,264 @@
+"""
+Closed-form solutions for the straight-shot limit of the cloth model.
+
+These are the textbook results the numerical integrator in :mod:`cueai.physics.ball`
+must reproduce, so they serve two purposes:
+
+1. Ground truth for the validation suite (``tests/test_validation.py``).
+2. A microsecond-cost baseline predictor that the learned residual model
+ corrects (:mod:`cueai.ml`), which is what makes the "physics-informed"
+ framing measurable rather than decorative.
+
+Derivation (ball of mass m, radius R, I = 2/5 mR², struck along +x)
+-------------------------------------------------------------------
+While the contact point slips, cloth friction acts opposite the slip velocity
+u with magnitude μ_s g, and the resulting torque drives u to zero at
+
+ |du/dt| = μ_s g (1 + mR²/I) = 3.5 μ_s g
+
+A tip offset of f·R above centre launches the ball with ω = 2.5 f v₀ / R, so
+u₀ = v₀ |1 - 2.5 f|:
+
+ t_slide = u₀ / (3.5 μ_s g)
+ v_roll = v₀ - s μ_s g t_slide (s = sign(1 - 2.5f))
+ d_slide = v₀ t_slide - s ½ μ_s g t_slide²
+ d_roll = v_roll² / (2 μ_r g)
+
+For a centre-ball hit (f = 0) this collapses to the familiar
+v_roll = 5/7 v₀ and d_slide = 12 v₀² / (49 μ_s g).
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+
+from cueai.physics.constants import BallParams, G, ShotParams, TableParams
+
+# Slip decays 3.5x faster than the centre of mass because the friction torque
+# also spins the ball up: 1 + mR²/I = 1 + 5/2.
+SLIP_DECAY_FACTOR = 3.5
+# Fraction of initial speed remaining once a centre-ball hit starts rolling.
+ROLLING_SPEED_RATIO = 5.0 / 7.0
+
+
+@dataclass(frozen=True)
+class StraightShot:
+ """Phase-resolved solution for a straight shot on an empty table."""
+
+ v_roll: float # speed at the sliding → rolling transition (m/s)
+ t_slide: float # duration of the sliding phase (s)
+ d_slide: float # distance covered while sliding (m)
+ d_roll: float # distance covered while rolling (m)
+
+ @property
+ def d_total(self) -> float:
+ return self.d_slide + self.d_roll
+
+
+def straight_shot(v0: float, table: TableParams, english_y: float = 0.0) -> StraightShot:
+ """Closed-form slide/roll decomposition for a straight shot."""
+ if v0 <= 0:
+ return StraightShot(0.0, 0.0, 0.0, 0.0)
+
+ slip0 = v0 * abs(1.0 - 2.5 * english_y)
+ sign = float(np.sign(1.0 - 2.5 * english_y))
+ t_slide = slip0 / (SLIP_DECAY_FACTOR * table.mu_slide * G)
+ v_roll = v0 - sign * table.mu_slide * G * t_slide
+ d_slide = v0 * t_slide - sign * 0.5 * table.mu_slide * G * t_slide**2
+ d_roll = v_roll**2 / (2.0 * table.mu_roll * G)
+ return StraightShot(
+ v_roll=float(v_roll),
+ t_slide=float(t_slide),
+ d_slide=float(d_slide),
+ d_roll=float(d_roll),
+ )
+
+
+def slide_distance(v0: float, mu_slide: float) -> float:
+ """Distance to the rolling transition for a centre-ball hit: 12v₀²/(49μg)."""
+ return 12.0 * v0**2 / (49.0 * mu_slide * G)
+
+
+def roll_distance(v_roll: float, mu_roll: float) -> float:
+ """Distance a rolling ball covers before stopping: v²/(2μ_r g)."""
+ return v_roll**2 / (2.0 * mu_roll * G)
+
+
+MAX_EVENTS = 32
+# Spatial resolution used to test a segment against the pocket mouths.
+POCKET_SCAN_STEP = 0.01
+
+
+@dataclass(frozen=True)
+class FreeBallOutcome:
+ """Where an unobstructed ball ends up, and what happened on the way."""
+
+ position: np.ndarray
+ potted: bool
+ cushions: int
+
+
+def _slide_rail_time(
+ pos: np.ndarray,
+ vel: np.ndarray,
+ accel: np.ndarray,
+ lo: np.ndarray,
+ hi: np.ndarray,
+ t_end: float,
+) -> tuple[float, int]:
+ """
+ First cushion crossing during a parabolic sliding segment.
+
+ Solves ½a t² + v t + (p - bound) = 0 per axis and boundary, keeping the
+ earliest root inside (0, t_end].
+ """
+ best_t, best_axis = t_end, -1
+ for axis in (0, 1):
+ for bound in (lo[axis], hi[axis]):
+ a, b, c = 0.5 * accel[axis], vel[axis], pos[axis] - bound
+ if abs(a) < 1e-14:
+ roots = [-c / b] if abs(b) > 1e-14 else []
+ else:
+ disc = b * b - 4 * a * c
+ if disc < 0:
+ continue
+ sqrt_disc = float(np.sqrt(disc))
+ roots = [(-b - sqrt_disc) / (2 * a), (-b + sqrt_disc) / (2 * a)]
+ for root in roots:
+ if 1e-9 < root <= best_t:
+ best_t, best_axis = float(root), axis
+ return best_t, best_axis
+
+
+def _ray_rail_distance(
+ pos: np.ndarray, direction: np.ndarray, lo: np.ndarray, hi: np.ndarray
+) -> tuple[float, int]:
+ """Distance to the first cushion along a straight rolling segment."""
+ best_distance, best_axis = np.inf, -1
+ for axis in (0, 1):
+ component = direction[axis]
+ if abs(component) < 1e-12:
+ continue
+ bound = hi[axis] if component > 0 else lo[axis]
+ distance = (bound - pos[axis]) / component
+ if 1e-9 < distance < best_distance:
+ best_distance, best_axis = float(distance), axis
+ return best_distance, best_axis
+
+
+def _first_pocket_hit(points: np.ndarray, table: TableParams) -> int:
+ """Index of the first sampled point inside a pocket mouth, or -1."""
+ if not table.pockets or len(points) == 0:
+ return -1
+ mouths = np.asarray(table.pockets, dtype=np.float64)
+ distances = np.linalg.norm(points[:, None, :] - mouths[None, :, :], axis=2)
+ inside = np.nonzero((distances < table.pocket_radius).any(axis=1))[0]
+ return int(inside[0]) if len(inside) else -1
+
+
+def _sample_count(length: float) -> int:
+ return int(np.clip(np.ceil(abs(length) / POCKET_SCAN_STEP), 2, 4096))
+
+
+def solve_free_ball(
+ pos: np.ndarray,
+ vel: np.ndarray,
+ spin_contact: np.ndarray,
+ table: TableParams,
+ radius: float = 0.028575,
+) -> FreeBallOutcome:
+ """
+ Closed-form trajectory of a single ball on an otherwise empty table.
+
+ Integrates nothing. The motion is a sequence of exactly solvable segments:
+
+ * **Sliding.** The slip velocity u decays along a fixed direction, so the
+ friction force is constant and the path is a parabola of known duration
+ |u| / (3.5 μ_s g).
+ * **Rolling.** A straight line of length v² / (2 μ_r g).
+ * **Cushion.** The normal velocity component reverses and loses energy while
+ the spin term u - v carries through, which is what makes a ball arrive at
+ the next rail sliding rather than rolling.
+
+ ``spin_contact`` is ω × (-Rẑ), the spin contribution to the contact-point
+ velocity, so the slip velocity is ``vel + spin_contact``.
+
+ Not modelled: ball-ball contact, rail friction and spin transfer, the
+ speed-dependent softening of the cushions, and cloth inhomogeneity.
+ """
+ pos = np.asarray(pos, dtype=np.float64).copy()
+ vel = np.asarray(vel, dtype=np.float64).copy()
+ spin = np.asarray(spin_contact, dtype=np.float64).copy()
+ lo = np.array([radius, radius])
+ hi = np.array([table.length - radius, table.width - radius])
+ cushions = 0
+
+ for _ in range(MAX_EVENTS):
+ slip = vel + spin
+ slip_speed = float(np.linalg.norm(slip))
+ speed = float(np.linalg.norm(vel))
+
+ if slip_speed > 1e-4:
+ accel = -table.mu_slide * G * (slip / slip_speed)
+ t_slide = slip_speed / (SLIP_DECAY_FACTOR * table.mu_slide * G)
+ t_hit, axis = _slide_rail_time(pos, vel, accel, lo, hi, t_slide)
+ t = min(t_slide, t_hit) if axis >= 0 else t_slide
+
+ steps = np.linspace(0.0, t, _sample_count(speed * t))[:, None]
+ path = pos + vel * steps + 0.5 * accel * steps**2
+ hit = _first_pocket_hit(path, table)
+ if hit >= 0:
+ return FreeBallOutcome(path[hit], True, cushions)
+
+ pos = pos + vel * t + 0.5 * accel * t**2
+ vel = vel + accel * t
+ spin = (slip - SLIP_DECAY_FACTOR * table.mu_slide * G * (slip / slip_speed) * t) - vel
+ if axis < 0 or t_hit >= t_slide:
+ spin = -vel # slip exhausted: the ball is now rolling
+ continue
+ else:
+ if speed < 1e-4:
+ return FreeBallOutcome(pos, False, cushions)
+ direction = vel / speed
+ stop_distance = roll_distance(speed, table.mu_roll)
+ rail_distance, axis = _ray_rail_distance(pos, direction, lo, hi)
+ travel = min(stop_distance, rail_distance)
+
+ steps = np.linspace(0.0, travel, _sample_count(travel))[:, None]
+ path = pos + direction * steps
+ hit = _first_pocket_hit(path, table)
+ if hit >= 0:
+ return FreeBallOutcome(path[hit], True, cushions)
+
+ pos = pos + direction * travel
+ if stop_distance <= rail_distance:
+ return FreeBallOutcome(pos, False, cushions)
+ speed = float(np.sqrt(max(speed**2 - 2 * table.mu_roll * G * travel, 0.0)))
+ vel = speed * direction
+ spin = -vel
+
+ # Cushion impulse: reverse and damp the normal velocity component. The
+ # spin term is untouched, so the ball leaves the rail sliding.
+ vel[axis] *= -table.e_cushion
+ pos[axis] = float(np.clip(pos[axis], lo[axis], hi[axis]))
+ cushions += 1
+
+ return FreeBallOutcome(pos, False, cushions)
+
+
+def predict_endpoint(
+ shot: ShotParams,
+ cue_pos: np.ndarray | tuple[float, float],
+ table: TableParams,
+ radius: float = 0.028575,
+) -> np.ndarray:
+ """Closed-form resting position of the cue ball for a given shot."""
+ omega = shot.initial_omega(BallParams(radius=radius))
+ # ω × (-Rẑ) = (-R ω_y, R ω_x)
+ spin_contact = np.array([-radius * omega[1], radius * omega[0]])
+ velocity = shot.speed * np.array([np.cos(shot.angle), np.sin(shot.angle)])
+ return solve_free_ball(
+ np.asarray(cue_pos, dtype=np.float64), velocity, spin_contact, table, radius
+ ).position
diff --git a/src/cueai/physics/ball.py b/src/cueai/physics/ball.py
index 72ca62b..5567de1 100644
--- a/src/cueai/physics/ball.py
+++ b/src/cueai/physics/ball.py
@@ -2,13 +2,14 @@
from __future__ import annotations
+import math
from dataclasses import dataclass, field
from enum import Enum, auto
-from typing import TYPE_CHECKING, Optional
+from typing import TYPE_CHECKING
import numpy as np
-from cueai.physics.constants import G, BallParams, TableParams
+from cueai.physics.constants import BallParams, G, TableParams
if TYPE_CHECKING:
from cueai.physics.rack import BallIdentity
@@ -31,9 +32,9 @@ class Ball:
params: BallParams = field(default_factory=BallParams)
pocketed: bool = False
number: int = 0
- identity: Optional["BallIdentity"] = None
+ identity: BallIdentity | None = None
- def copy(self) -> "Ball":
+ def copy(self) -> Ball:
return Ball(
id=self.id,
pos=self.pos.copy(),
@@ -58,11 +59,16 @@ def I(self) -> float:
return self.params.inertia
def slip_velocity(self) -> np.ndarray:
- """Velocity of cloth contact point: u = v + ω × (-Rẑ)."""
+ """
+ Velocity of the cloth contact point: u = v + ω × (-Rẑ).
+
+ Expanding the cross product gives u = (vₓ - Rω_y, v_y + Rω_x), so pure
+ rolling (u = 0) corresponds to ω_y = vₓ/R and ω_x = -v_y/R.
+ """
return np.array(
[
- self.vel[0] + self.omega[1] * self.R,
- self.vel[1] - self.omega[0] * self.R,
+ self.vel[0] - self.omega[1] * self.R,
+ self.vel[1] + self.omega[0] * self.R,
],
dtype=np.float64,
)
@@ -77,7 +83,9 @@ def motion_state(self, table: TableParams, eps: float = 1e-3) -> MotionState:
u_mag = float(np.linalg.norm(u))
v_mag = self.speed()
wz = abs(self.omega[2])
- slip_eps = max(0.05, 0.05 * max(v_mag, 0.1))
+ # Tolerance must exceed the per-step slip decrement (≈3.5 μ g dt) to avoid
+ # chattering between SLIDING and ROLLING near the transition.
+ slip_eps = max(0.01, 0.01 * v_mag)
if v_mag < eps and u_mag < eps and wz < eps:
return MotionState.STATIONARY
if u_mag < slip_eps and v_mag >= eps:
@@ -87,6 +95,22 @@ def motion_state(self, table: TableParams, eps: float = 1e-3) -> MotionState:
return MotionState.SLIDING
+def decay_spin(omega_z: float, table: TableParams, radius: float, dt: float) -> float:
+ """
+ Spin-down about the vertical axis, clamped at zero.
+
+ The decrement per step is constant, so subtracting it unconditionally steps
+ straight past zero and flips the sign; the ball then chatters between two
+ small values and never reaches rest, and `Simulator.run` gives up at
+ `max_time` instead of stopping when the table is still. Clamping is not a
+ tolerance hack: friction removes spin, it cannot reverse it.
+ """
+ step = 2.5 * table.mu_spin * G / radius * dt
+ if abs(omega_z) <= step:
+ return 0.0
+ return omega_z - math.copysign(step, omega_z)
+
+
def local_mu_slide(table: TableParams, pos: np.ndarray) -> float:
"""Spatially varying sliding friction (table imperfections via smooth noise)."""
amp = table.friction_noise_amp
@@ -108,7 +132,6 @@ def integrate_ball(ball: Ball, table: TableParams, dt: float) -> Ball:
R, m, I = b.R, b.m, b.I
mu_s = local_mu_slide(table, b.pos)
mu_r = table.mu_roll
- mu_sp = table.mu_spin
if state is MotionState.STATIONARY:
b.vel[:] = 0
@@ -116,11 +139,7 @@ def integrate_ball(ball: Ball, table: TableParams, dt: float) -> Ball:
return b
if state is MotionState.SPINNING:
- decay = 2.5 * mu_sp * G / R
- sgn = np.sign(b.omega[2]) if b.omega[2] != 0 else 0.0
- b.omega[2] -= decay * sgn * dt
- if abs(b.omega[2]) < 1e-4:
- b.omega[2] = 0.0
+ b.omega[2] = decay_spin(float(b.omega[2]), table, R, dt)
return b
if state is MotionState.SLIDING:
@@ -129,14 +148,14 @@ def integrate_ball(ball: Ball, table: TableParams, dt: float) -> Ball:
if u_mag < 1e-9:
return b
u_hat = u / u_mag
- # Curving force: slip has a component from sidespin → path curves
+ # Friction opposes the contact-point slip, not the ball centre velocity
a = -mu_s * G * u_hat
- alpha = np.array([-R * m * a[1] / I, R * m * a[0] / I, 0.0])
- if abs(b.omega[2]) > 1e-6:
- alpha[2] = -np.sign(b.omega[2]) * 2.5 * mu_sp * G / R
+ # τ = (-Rẑ) × F ⇒ α = (R m a_y / I, -R m a_x / I, 0)
+ alpha = np.array([R * m * a[1] / I, -R * m * a[0] / I, 0.0])
b.vel = b.vel + a * dt
b.omega = b.omega + alpha * dt
+ b.omega[2] = decay_spin(float(b.omega[2]), table, R, dt)
u2 = b.slip_velocity()
if float(np.linalg.norm(u2)) < max(1e-3, 0.01 * b.speed()):
b.omega[0] = -b.vel[1] / R
@@ -150,12 +169,9 @@ def integrate_ball(ball: Ball, table: TableParams, dt: float) -> Ball:
b.omega[:2] = 0
return b
v_hat = b.vel / v_mag
- # Mild curve while rolling if residual ωz (english hold)
+ # Rolling resistance only. Vertical-axis spin exerts no lateral force on a
+ # rolling rigid sphere; it acts through cushion and ball-ball throw instead.
a = -mu_r * G * v_hat
- if abs(b.omega[2]) > 0.5:
- # small lateral force from residual spin-on-cloth coupling
- lateral = np.array([-v_hat[1], v_hat[0]]) * (0.002 * b.omega[2])
- a = a + lateral
b.vel = b.vel + a * dt
if b.speed() < 2e-2:
b.vel[:] = 0
@@ -163,7 +179,6 @@ def integrate_ball(ball: Ball, table: TableParams, dt: float) -> Ball:
return b
b.omega[0] = -b.vel[1] / R
b.omega[1] = b.vel[0] / R
- if abs(b.omega[2]) > 1e-6:
- b.omega[2] -= np.sign(b.omega[2]) * 2.5 * mu_sp * G / R * dt
+ b.omega[2] = decay_spin(float(b.omega[2]), table, R, dt)
b.pos = b.pos + b.vel * dt
return b
diff --git a/src/cueai/physics/collisions.py b/src/cueai/physics/collisions.py
index 10db534..4fce6c5 100644
--- a/src/cueai/physics/collisions.py
+++ b/src/cueai/physics/collisions.py
@@ -7,6 +7,19 @@
from cueai.physics.ball import Ball
from cueai.physics.constants import TableParams
+# Two balls count as touching once their surfaces are inside this band.
+#
+# The value is a constant rather than a literal because it has a correctness
+# requirement, not just a tuning one: it must sit far above the floating-point
+# noise in a position (~1e-16 m) and far below anything physical, and it must
+# not coincide with the gap balls are racked at. It used to be 1e-4, exactly
+# the racking clearance, and the consequence was that only sixteen of a rack's
+# thirty contacts registered — which sixteen being decided by whether hypot()
+# happened to round up or down. A break then propagated through a contact graph
+# with holes in it, so balls in the middle of the rack came out of a full-power
+# break having never moved.
+CONTACT_BAND = 1e-5
+
def ball_ball_friction(v_rel: float, table: TableParams) -> float:
"""
@@ -25,24 +38,31 @@ def resolve_ball_ball(a: Ball, b: Ball, table: TableParams) -> tuple[Ball, Ball]
Includes relative-velocity-dependent friction and contact-point ω coupling.
"""
- a, b = a.copy(), b.copy()
if a.pocketed or b.pocketed:
return a, b
delta = b.pos - a.pos
- dist = float(np.linalg.norm(delta))
+ dist = float(np.hypot(delta[0], delta[1]))
min_dist = a.R + b.R
- if dist < 1e-12:
- return a, b
- if dist > min_dist + 2e-4:
+ if dist < 1e-12 or dist > min_dist + CONTACT_BAND:
return a, b
n = delta / dist
overlap = min_dist - dist
+ # Angular velocity contributes nothing along n, because ω × (R n) ⊥ n, so the
+ # approach speed can be tested before doing any of the impulse work. Balls
+ # resting in contact — most pairs in a packed rack, every step — stop here.
+ v_n = float((a.vel - b.vel) @ n)
+ if v_n <= 1e-6 and overlap <= 0:
+ return a, b
+
+ a, b = a.copy(), b.copy()
if overlap > 0:
# Mass-proportional separation (equal mass → half/half)
a.pos = a.pos - 0.5 * overlap * n
b.pos = b.pos + 0.5 * overlap * n
+ if v_n <= 1e-6:
+ return a, b # separating or resting: positional correction only
ra = a.R * np.array([n[0], n[1], 0.0])
rb = -b.R * np.array([n[0], n[1], 0.0])
@@ -52,11 +72,8 @@ def resolve_ball_ball(a: Ball, b: Ball, table: TableParams) -> tuple[Ball, Ball]
vb3 = np.array([b.vel[0], b.vel[1], 0.0])
v_rel = (va3 - vb3) + np.cross(wa, ra) - np.cross(wb, rb)
- v_n = float(np.dot(v_rel[:2], n))
- # n points a→b; approaching when (va−vb)·n > 0
- if v_n <= 1e-6:
- return a, b # separating or resting
-
+ # v_n was already established above from the linear velocities; the spin terms
+ # in v_rel are tangential and only matter for the friction impulse below.
e = table.e_ball
# Slightly softer at high speed (energy loss grows)
speed_n = abs(v_n)
@@ -183,34 +200,57 @@ def check_pocket(ball: Ball, table: TableParams) -> Ball:
def resolve_all_ball_collisions(
- balls: list[Ball], table: TableParams, passes: int = 8
+ balls: list[Ball], table: TableParams, passes: int = 24
) -> list[Ball]:
"""
Multi-pass pairwise resolution so cluster breaks / simultaneous contacts
propagate (critical for a packed rack).
+
+ The count is a safety cap, not the usual cost: the loop exits as soon as a
+ pass changes nothing, which for a table that is merely resting is the first
+ one. It has to exceed the depth of the contact chain the impulse travels
+ along, five rows for a full rack, and `web/js/physics.js` must use the same
+ number or the two implementations diverge on exactly the shots that matter.
"""
balls = [b.copy() for b in balls]
n = len(balls)
for _ in range(passes):
any_hit = False
+ # Candidate search is vectorised: with 16 balls this runs every time step,
+ # so the pairwise distances are one numpy call rather than 120 Python ones.
+ active = [i for i in range(n) if not balls[i].pocketed]
+ if len(active) < 2:
+ break
+ positions = np.array([balls[i].pos for i in active])
+ radii = np.array([balls[i].R for i in active])
+ gaps = (
+ np.linalg.norm(positions[:, None, :] - positions[None, :, :], axis=2)
+ - radii[:, None]
+ - radii[None, :]
+ )
+ rows, cols = np.nonzero(np.triu(gaps <= CONTACT_BAND, k=1))
# Process deepest overlaps first for stability
- pairs: list[tuple[float, int, int]] = []
- for i in range(n):
- if balls[i].pocketed:
- continue
- for j in range(i + 1, n):
- if balls[j].pocketed:
- continue
- d = float(np.linalg.norm(balls[i].pos - balls[j].pos))
- min_d = balls[i].R + balls[j].R
- if d <= min_d + 1e-4:
- pairs.append((d - min_d, i, j))
- pairs.sort() # most overlapped first
+ pairs = sorted(
+ (float(gaps[r, c]), active[r], active[c]) for r, c in zip(rows, cols)
+ )
for _, i, j in pairs:
bi, bj = resolve_ball_ball(balls[i], balls[j], table)
- if not np.allclose(bi.vel, balls[i].vel) or not np.allclose(bj.vel, balls[j].vel):
+ # Exact comparison, not ``allclose``: its 1e-8 absolute tolerance
+ # would end the sweep a pass early on nanometre-scale corrections,
+ # and since the JavaScript port reports what it did rather than
+ # inferring it, a tolerance here is a divergence between the two.
+ if (
+ not np.array_equal(bi.vel, balls[i].vel)
+ or not np.array_equal(bj.vel, balls[j].vel)
+ or not np.array_equal(bi.pos, balls[i].pos)
+ or not np.array_equal(bj.pos, balls[j].pos)
+ ):
any_hit = True
balls[i], balls[j] = bi, bj
- if not any_hit and not pairs:
+ # A pass that changed nothing leaves the next one the same problem, so
+ # this is convergence rather than a budget. It matters now that a
+ # racked table reports thirty resting contacts every step instead of
+ # sixteen: without it every step below would pay for all the passes.
+ if not any_hit:
break
return balls
diff --git a/src/cueai/physics/constants.py b/src/cueai/physics/constants.py
index 9470f33..9b33738 100644
--- a/src/cueai/physics/constants.py
+++ b/src/cueai/physics/constants.py
@@ -6,7 +6,6 @@
import numpy as np
-
G = 9.81 # m/s^2
@@ -65,17 +64,32 @@ class ShotParams:
speed: float # tip speed imparted → cue-ball linear speed (m/s)
angle: float # launch angle radians, 0 = +x
- english_x: float = 0 # sidespin tip offset as fraction of R (-1..1)
- english_y: float = 0 # topspin / backspin tip offset (-1..1)
+ english_x: float = 0 # sidespin tip offset as fraction of R (+ = right english)
+ english_y: float = 0 # topspin (+) / backspin (−) tip offset (-1..1)
cue_elevation: float = 0.0 # radians (jump / massé lite)
def initial_omega(self, ball: BallParams) -> np.ndarray:
- """Map tip offset to angular velocity (simplified cue model)."""
+ """
+ Map cue tip offset to initial angular velocity.
+
+ A horizontal impulse J applied a distance d = f·R off centre gives
+ Δv = J/m and Δω = J·d/I, so with I = (2/5)mR²:
+
+ ω = 2.5 · f · v / R
+
+ Hence f = 0.4 (tip 2R/5 above centre, i.e. 7R/5 above the cloth)
+ launches the ball already rolling — the standard "natural roll" result —
+ and |f| > 0.5 is past the practical miscue limit.
+
+ - Top/backspin acts about the horizontal axis perpendicular to travel,
+ which for a shot along θ is (-sin θ, cos θ, 0) — the same axis as
+ natural roll, so english_y = +1 is pure follow.
+ - Sidespin acts about the vertical axis; positive english_x (right
+ english) is clockwise seen from above, hence negative ω_z.
+ """
R = ball.radius
- # Tip contact ≈ tangential impulse → ω ≈ (v_tip × r) / (k R)
- # Use Alciatore-style: ω_y ~ -english_x * v / R, ω_x ~ english_y * v / R
v = self.speed
- wx = self.english_y * v / R
- wy = -self.english_x * v / R
- wz = 0.15 * self.english_x * v / R # slight vertical-axis from sidespin
- return np.array([wx, wy, wz], dtype=np.float64)
+ roll_axis = np.array([-np.sin(self.angle), np.cos(self.angle), 0.0])
+ omega = 2.5 * self.english_y * (v / R) * roll_axis
+ omega[2] -= 2.5 * self.english_x * v / R
+ return omega
diff --git a/src/cueai/physics/rack.py b/src/cueai/physics/rack.py
index a32283b..901d8d1 100644
--- a/src/cueai/physics/rack.py
+++ b/src/cueai/physics/rack.py
@@ -69,7 +69,10 @@ def triangle_positions(apex: np.ndarray, R: float) -> list[np.ndarray]:
15 tightly packed positions: row 0 has 1 ball (apex), … row 4 has 5.
Oriented toward the head of the table (−x).
"""
- gap = 2.0 * R + 1e-4 # tiny clearance so they aren't already overlapping
+ # Racked balls touch. The clearance that used to be here was the same 1e-4
+ # as the collision solver's contact band, which put every contact in the
+ # rack exactly on the threshold that decides whether it exists.
+ gap = 2.0 * R
positions: list[np.ndarray] = []
for row in range(5):
for col in range(row + 1):
@@ -142,7 +145,7 @@ def make_full_rack(
)
balls.append(cue)
- for idx, (num, pos) in enumerate(zip(numbers, spots)):
+ for num, pos in zip(numbers, spots):
balls.append(
Ball(
id=num, # use ball number as id for clarity
diff --git a/src/cueai/physics/simulator.py b/src/cueai/physics/simulator.py
index 1280d29..fb7344a 100644
--- a/src/cueai/physics/simulator.py
+++ b/src/cueai/physics/simulator.py
@@ -27,6 +27,7 @@ class SimResult:
endpoints: dict[int, np.ndarray]
ball_meta: dict[int, dict] = field(default_factory=dict)
collision_events: int = 0
+ cushion_events: int = 0
@dataclass
@@ -103,6 +104,7 @@ def run(self, balls: list[Ball], record_every: int = 4) -> SimResult:
vels: dict[int, list[np.ndarray]] = {b.id: [] for b in balls}
omgs: dict[int, list[np.ndarray]] = {b.id: [] for b in balls}
collisions = 0
+ cushions = 0
rest_steps = 0
meta = {
@@ -121,7 +123,10 @@ def run(self, balls: list[Ball], record_every: int = 4) -> SimResult:
if b.pocketed:
continue
balls[i] = integrate_ball(b, self.table, self.dt)
+ pre_rail_vel = balls[i].vel.copy()
balls[i] = resolve_cushion(balls[i], self.table)
+ if not np.allclose(pre_rail_vel, balls[i].vel, atol=1e-6):
+ cushions += 1
balls[i] = check_pocket(balls[i], self.table)
if balls[i].speed() > 1e-4 or float(np.linalg.norm(balls[i].omega)) > 1e-3:
moving = True
@@ -197,6 +202,7 @@ def run(self, balls: list[Ball], record_every: int = 4) -> SimResult:
endpoints=endpoints,
ball_meta=meta,
collision_events=collisions,
+ cushion_events=cushions,
)
def simulate_shot(
diff --git a/src/cueai/ui/app.py b/src/cueai/ui/app.py
index a70fd3e..2280bcb 100644
--- a/src/cueai/ui/app.py
+++ b/src/cueai/ui/app.py
@@ -7,20 +7,20 @@
from pathlib import Path
import numpy as np
-from PyQt6.QtCore import QPointF, Qt, QTimer, QRectF
-from PyQt6.QtGui import QColor, QPainter, QPen, QBrush, QFont
+from PyQt6.QtCore import QPointF, QRectF, Qt, QTimer
+from PyQt6.QtGui import QBrush, QColor, QFont, QPainter, QPen
from PyQt6.QtWidgets import (
QApplication,
+ QCheckBox,
QDoubleSpinBox,
QFormLayout,
QHBoxLayout,
QLabel,
QMainWindow,
QPushButton,
- QCheckBox,
+ QSlider,
QVBoxLayout,
QWidget,
- QSlider,
)
from cueai.ml.infer import TrajectoryPredictor
@@ -339,7 +339,7 @@ def __init__(self):
self.noise.setRange(0, 0.08)
self.noise.setValue(0.025)
self.noise.setSingleStep(0.005)
- self.use_ml = QCheckBox("ML residual fusion")
+ self.use_ml = QCheckBox("Compare fast prediction")
self.use_ml.setChecked(True)
self.show_paths = QCheckBox("Show trails")
self.show_paths.setChecked(True)
@@ -408,9 +408,7 @@ def shoot(self) -> None:
mu_slide=self.mu.value(),
friction_noise_amp=self.noise.value(),
)
- # Continue from current table state (multi-shot play)
- live = [b.copy() for b in self.canvas.balls if not b.pocketed or b.number == 0]
- # keep pocketed markers
+ # Continue from the current table state, pocketed markers included
all_balls = [b.copy() for b in self.canvas.balls]
result = self.predictor.predict(
shot,
@@ -432,12 +430,20 @@ def _freeze() -> None:
sunk = [k for k, v in result["pocketed"].items() if v and k != "0"]
cue_scratch = bool(result["pocketed"].get("0"))
- ml = "ON" if result["ml_loaded"] and self.use_ml.isChecked() else "OFF"
- self.status.setText(
- f"ML {ml} · collisions≈{result.get('collisions', 0)} · "
+ timing = result.get("timing_ms", {})
+ gap = result.get("agreement", {}).get("cue_gap_m")
+ summary = (
+ f"simulated in {timing.get('simulator', 0.0) / 1000:.1f} s · "
+ f"collisions≈{result.get('collisions', 0)} · cushions≈{result.get('cushions', 0)} · "
f"pocketed {len(sunk)} {sunk[:6]}{'…' if len(sunk) > 6 else ''}"
+ (" · SCRATCH" if cue_scratch else "")
)
+ if self.use_ml.isChecked() and gap is not None:
+ summary += (
+ f"\nfast prediction took {timing.get('fast', 0.0):.2f} ms and landed "
+ f"{gap * 1000:.0f} mm from the simulated cue ball"
+ )
+ self.status.setText(summary)
def main() -> None:
diff --git a/tests/test_api.py b/tests/test_api.py
new file mode 100644
index 0000000..30d324b
--- /dev/null
+++ b/tests/test_api.py
@@ -0,0 +1,80 @@
+"""Contract tests for the FastAPI surface, driven through the ASGI app directly."""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+from fastapi.testclient import TestClient
+
+from cueai.api.main import app
+
+SOFT_SHOT = {"speed": 1.5, "angle_deg": 5.0, "full_rack": False, "cue_x": 0.6, "cue_y": 0.6}
+
+
+@pytest.fixture(scope="module")
+def client() -> TestClient:
+ return TestClient(app)
+
+
+def test_health_reports_model_state(client: TestClient) -> None:
+ body = client.get("/health").json()
+ assert body["status"] == "ok"
+ assert body["backend"] in {"none", "torch", "onnx"}
+ assert isinstance(body["ml_loaded"], bool)
+
+
+def test_table_geometry_is_a_nine_foot_table(client: TestClient) -> None:
+ body = client.get("/table").json()
+ assert body["length"] == pytest.approx(2.54)
+ assert body["width"] == pytest.approx(1.27)
+ assert len(body["pockets"]) == 6
+
+
+def test_rack_is_a_legal_eight_ball_rack(client: TestClient) -> None:
+ balls = client.get("/rack", params={"seed": 7}).json()["balls"]
+ assert len(balls) == 16
+ assert sorted(b["number"] for b in balls) == list(range(16))
+ eight = next(b for b in balls if b["number"] == 8)
+ assert eight["suit"] == "eight"
+
+
+def test_fast_prediction_stays_on_the_table(client: TestClient) -> None:
+ body = client.post("/predict/fast", json=SOFT_SHOT).json()
+ cue = body["endpoints"]["cue"]
+ assert 0.0 <= cue[0] <= 2.54
+ assert 0.0 <= cue[1] <= 1.27
+ assert body["timing_ms"] < 50
+
+
+def test_fast_prediction_omits_the_object_ball_when_none_was_given(client: TestClient) -> None:
+ """Reporting a position for a ball the caller never placed would be noise."""
+ body = client.post("/predict/fast", json=SOFT_SHOT).json()
+ assert body["endpoints"]["object"] is None
+
+ with_object = client.post("/predict/fast", json={**SOFT_SHOT, "obj_x": 1.4, "obj_y": 0.7})
+ assert len(with_object.json()["endpoints"]["object"]) == 2
+
+
+def test_fast_prediction_is_reproducible(client: TestClient) -> None:
+ first = client.post("/predict/fast", json=SOFT_SHOT).json()
+ second = client.post("/predict/fast", json=SOFT_SHOT).json()
+ assert first["endpoints"] == second["endpoints"]
+
+
+def test_simulation_returns_trajectories_and_timings(client: TestClient) -> None:
+ body = client.post("/predict", json={**SOFT_SHOT, "obj_x": 1.4, "obj_y": 0.7}).json()
+ cue_path = np.asarray(body["trajectory"]["0"])
+ assert cue_path.shape[1] == 2
+ assert len(cue_path) > 10
+ assert body["timing_ms"]["simulator"] > body["timing_ms"]["fast"]
+
+
+def test_rejects_out_of_range_shots(client: TestClient) -> None:
+ assert client.post("/predict/fast", json={**SOFT_SHOT, "speed": 99.0}).status_code == 422
+ assert client.post("/predict/fast", json={**SOFT_SHOT, "english_x": 3.0}).status_code == 422
+
+
+def test_rejects_tip_offsets_past_the_miscue_limit(client: TestClient) -> None:
+ """The model never saw offsets beyond the miscue limit, so the API must not accept them."""
+ assert client.post("/predict/fast", json={**SOFT_SHOT, "english_y": 0.9}).status_code == 422
+ assert client.post("/predict/fast", json={**SOFT_SHOT, "english_y": 0.5}).status_code == 200
diff --git a/tests/test_features.py b/tests/test_features.py
new file mode 100644
index 0000000..df808a4
--- /dev/null
+++ b/tests/test_features.py
@@ -0,0 +1,100 @@
+"""
+Guards against train/serve skew.
+
+A residual model is only as good as the assumption that features are built the
+same way at training time and at request time. These tests pin that down, since
+a silent mismatch would show up as a quietly worse model rather than an error.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from cueai.ml import features as feat
+from cueai.ml.dataset import simulate_sample
+from cueai.physics.constants import ShotParams, TableParams
+
+
+def test_feature_names_match_vector_length() -> None:
+ vector = feat.build_features(
+ ShotParams(speed=2.0, angle=0.3, english_x=0.2, english_y=-0.1),
+ (0.6, 0.5),
+ (1.4, 0.7),
+ TableParams(),
+ )
+ assert vector.shape == (len(feat.FEATURE_NAMES),)
+ assert np.isfinite(vector).all()
+
+
+def test_serving_path_matches_training_path() -> None:
+ """
+ The dataset stores raw parameters and rebuilds features in bulk; serving
+ builds them one shot at a time. The two must agree exactly.
+ """
+ rows = [simulate_sample(index, seed=11) for index in range(4)]
+ df = pd.DataFrame(rows)
+ bulk = feat.build_feature_frame(df)
+
+ for position, record in enumerate(rows):
+ shot = ShotParams(
+ speed=record["speed"],
+ angle=record["angle"],
+ english_x=record["english_x"],
+ english_y=record["english_y"],
+ cue_elevation=record["cue_elevation"],
+ )
+ table = TableParams(
+ mu_slide=record["mu_slide"],
+ mu_roll=record["mu_roll"],
+ mu_spin=record["mu_spin"],
+ e_cushion=record["e_cushion"],
+ friction_noise_amp=record["friction_noise_amp"],
+ )
+ single = feat.build_features(
+ shot, (record["cue_x"], record["cue_y"]), (record["obj_x"], record["obj_y"]), table
+ )
+ np.testing.assert_allclose(single, bulk[position], rtol=0, atol=0)
+
+
+def test_baseline_features_agree_with_the_stored_baseline() -> None:
+ """The closed-form endpoint in the features is the one used as the residual base."""
+ record = simulate_sample(3, seed=11)
+ vector = feat.build_features(
+ ShotParams(
+ speed=record["speed"],
+ angle=record["angle"],
+ english_x=record["english_x"],
+ english_y=record["english_y"],
+ cue_elevation=record["cue_elevation"],
+ ),
+ (record["cue_x"], record["cue_y"]),
+ (record["obj_x"], record["obj_y"]),
+ TableParams(
+ mu_slide=record["mu_slide"],
+ mu_roll=record["mu_roll"],
+ mu_spin=record["mu_spin"],
+ e_cushion=record["e_cushion"],
+ friction_noise_amp=record["friction_noise_amp"],
+ ),
+ )
+ named = dict(zip(feat.FEATURE_NAMES, vector))
+ assert named["base_cue_x"] == pytest.approx(record["baseline_cue_end_x"])
+ assert named["base_cue_y"] == pytest.approx(record["baseline_cue_end_y"])
+
+
+def test_contact_geometry_flags_a_straight_on_shot() -> None:
+ table = TableParams()
+ straight = feat.build_features(
+ ShotParams(speed=2.0, angle=0.0), (0.6, 0.635), (1.2, 0.635), table
+ )
+ named = dict(zip(feat.FEATURE_NAMES, straight))
+ assert named["will_contact"] == 1.0
+ assert named["contact_perp"] == pytest.approx(0.0, abs=1e-9)
+ assert named["contact_along"] == pytest.approx(0.6)
+
+ wide = feat.build_features(
+ ShotParams(speed=2.0, angle=0.6), (0.6, 0.635), (1.2, 0.635), table
+ )
+ assert dict(zip(feat.FEATURE_NAMES, wide))["will_contact"] == 0.0
diff --git a/tests/test_metrics.py b/tests/test_metrics.py
new file mode 100644
index 0000000..a216de2
--- /dev/null
+++ b/tests/test_metrics.py
@@ -0,0 +1,54 @@
+"""
+Tests for the reported evaluation metrics.
+
+These functions produce the numbers the README quotes, so a silent error here
+would misstate the results rather than fail a build. The gate direction in
+particular is easy to invert without anything looking wrong.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from cueai.ml.train import endpoint_errors, risk_coverage
+
+
+def test_endpoint_errors_measures_distance_not_axes() -> None:
+ """A 3-4-5 offset on both balls is a 5 mm error, not a 3.5 mm mean of axes."""
+ truth = np.zeros((1, 4))
+ pred = np.array([[0.003, 0.004, 0.003, 0.004]])
+ errors = endpoint_errors(pred, truth)
+ assert errors["euclidean_mm"] == 5.0
+ assert errors["cue_mm"] == 5.0
+ assert errors["obj_mm"] == 5.0
+
+
+def test_risk_coverage_gates_on_the_expected_cushion_count() -> None:
+ """
+ Tighter gates must answer fewer shots and, when the gate carries signal,
+ answer them more accurately.
+ """
+ expected = np.repeat([0, 1, 2, 3], 10).astype(float)
+ truth = np.zeros((40, 4))
+ # Error grows with the expected cushion count, which is the pattern the
+ # real gate exploits.
+ pred = np.tile((expected * 0.1)[:, None], (1, 4))
+
+ rows = risk_coverage({"model": pred}, truth, expected)
+ labels = [row["expected_cushions_at_most"] for row in rows]
+ assert labels == ["0", "1", "2", "3", "all"]
+
+ coverage = [row["coverage_pct"] for row in rows]
+ assert coverage == sorted(coverage)
+ assert coverage[-1] == 100.0
+
+ errors = [row["model"] for row in rows]
+ assert errors == sorted(errors)
+ assert errors[0] == 0.0
+
+
+def test_risk_coverage_skips_thresholds_with_no_shots() -> None:
+ expected = np.full(5, 3.0)
+ truth = np.zeros((5, 4))
+ rows = risk_coverage({"model": truth}, truth, expected)
+ assert [row["expected_cushions_at_most"] for row in rows] == ["3", "all"]
diff --git a/tests/test_physics.py b/tests/test_physics.py
index e58fe7c..58c95a9 100644
--- a/tests/test_physics.py
+++ b/tests/test_physics.py
@@ -4,7 +4,7 @@
import numpy as np
-from cueai.physics.ball import Ball, integrate_ball, MotionState
+from cueai.physics.ball import Ball, MotionState, integrate_ball
from cueai.physics.constants import BallParams, ShotParams, TableParams
from cueai.physics.rack import make_full_rack
from cueai.physics.simulator import Simulator, shot_feature_vector
@@ -73,4 +73,8 @@ def test_collision_transfers_momentum():
obj_pos=(1.2, 0.635),
full_rack=False,
)
- assert result.endpoints[1][0] > 1.2
+ # The object ball starts at rest, so any speed at all came from the cue ball.
+ obj_speed = np.linalg.norm(result.velocities[1], axis=1)
+ cue_speed = np.linalg.norm(result.velocities[0], axis=1)
+ assert obj_speed.max() > 0.5 * cue_speed.max()
+ assert np.abs(np.diff(result.trajectories[1][:, 0])).sum() > 0.5
diff --git a/tests/test_validation.py b/tests/test_validation.py
new file mode 100644
index 0000000..481ccb6
--- /dev/null
+++ b/tests/test_validation.py
@@ -0,0 +1,516 @@
+"""
+Validation of the numerical simulator against closed-form physics.
+
+Every test here compares the integrator to an independently derived analytic
+result or to a conservation law, so a failure means the physics is wrong rather
+than merely different. The tolerances are the numbers quoted in
+``docs/VALIDATION.md``.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+from cueai.physics import analytic
+from cueai.physics.ball import Ball, MotionState, integrate_ball
+from cueai.physics.collisions import (
+ CONTACT_BAND,
+ resolve_all_ball_collisions,
+ resolve_ball_ball,
+ resolve_cushion,
+)
+from cueai.physics.constants import BallParams, G, ShotParams, TableParams
+from cueai.physics.rack import make_full_rack
+from cueai.physics.simulator import Simulator
+
+DT = 1e-4
+SMOOTH_CLOTH = {"friction_noise_amp": 0.0}
+
+
+def _launch(v0: float, english_y: float = 0.0, english_x: float = 0.0) -> Ball:
+ params = BallParams()
+ shot = ShotParams(speed=v0, angle=0.0, english_x=english_x, english_y=english_y)
+ return Ball(
+ id=0,
+ pos=np.array([0.0, 0.635]),
+ vel=np.array([v0, 0.0]),
+ omega=shot.initial_omega(params),
+ params=params,
+ )
+
+
+def _roll_out(
+ ball: Ball,
+ table: TableParams,
+ dt: float = DT,
+ max_time: float = 60.0,
+ stop_at_transition: bool = False,
+):
+ """Integrate an unobstructed ball to rest, capturing the rolling transition."""
+ transition: tuple[float, float] | None = None
+ for _ in range(int(max_time / dt)):
+ ball = integrate_ball(ball, table, dt)
+ if transition is None and ball.motion_state(table) is MotionState.ROLLING:
+ transition = (ball.speed(), float(ball.pos[0]))
+ if stop_at_transition:
+ break
+ if ball.speed() == 0.0:
+ break
+ return ball, transition
+
+
+@pytest.mark.parametrize("v0", [1.0, 2.0, 4.0])
+def test_rolling_speed_is_five_sevenths_of_launch_speed(v0: float) -> None:
+ """Classic result: a centre-ball hit starts rolling at 5/7 v₀."""
+ table = TableParams(**SMOOTH_CLOTH)
+ _, transition = _roll_out(_launch(v0), table, stop_at_transition=True)
+ assert transition is not None
+ v_roll, _ = transition
+ assert v_roll == pytest.approx(analytic.ROLLING_SPEED_RATIO * v0, rel=0.01)
+
+
+@pytest.mark.parametrize("v0", [1.0, 2.0, 4.0])
+def test_slide_distance_matches_closed_form(v0: float) -> None:
+ """Sliding phase length must equal 12v₀²/(49 μ_s g)."""
+ table = TableParams(**SMOOTH_CLOTH)
+ _, transition = _roll_out(_launch(v0), table, stop_at_transition=True)
+ assert transition is not None
+ _, x_roll = transition
+ assert x_roll == pytest.approx(analytic.slide_distance(v0, table.mu_slide), rel=0.02)
+
+
+@pytest.mark.parametrize("v0", [1.0, 2.0, 3.0])
+def test_total_stopping_distance_matches_closed_form(v0: float) -> None:
+ """Slide + roll distance on an unbounded table, within 1%."""
+ # A brisk cloth (mu_roll=0.03) keeps the integration short enough for CI.
+ table = TableParams(mu_roll=0.03, **SMOOTH_CLOTH)
+ ball, _ = _roll_out(_launch(v0), table)
+ expected = analytic.straight_shot(v0, table).d_total
+ assert float(ball.pos[0]) == pytest.approx(expected, rel=0.01)
+
+
+def test_natural_roll_tip_offset_produces_no_sliding_phase() -> None:
+ """A tip 2R/5 above centre launches the ball already rolling."""
+ table = TableParams(mu_roll=0.03, **SMOOTH_CLOTH)
+ ball = _launch(2.0, english_y=0.4)
+ assert ball.motion_state(table) is MotionState.ROLLING
+ assert float(np.linalg.norm(ball.slip_velocity())) < 1e-9
+ stopped, _ = _roll_out(ball, table)
+ expected = analytic.roll_distance(2.0, table.mu_roll)
+ assert float(stopped.pos[0]) == pytest.approx(expected, rel=0.01)
+ assert float(stopped.pos[1]) == pytest.approx(0.635, abs=1e-9)
+
+
+def test_spin_decays_at_analytic_rate() -> None:
+ """Vertical-axis spin decays linearly at 5 μ_sp g / 2R."""
+ table = TableParams(**SMOOTH_CLOTH)
+ params = BallParams()
+ omega0 = 60.0
+ ball = Ball(
+ id=0,
+ pos=np.array([1.0, 0.6]),
+ vel=np.zeros(2),
+ omega=np.array([0.0, 0.0, omega0]),
+ params=params,
+ )
+ steps = 500
+ for _ in range(steps):
+ ball = integrate_ball(ball, table, DT)
+ expected = omega0 - 2.5 * table.mu_spin * G / params.radius * steps * DT
+ assert float(ball.omega[2]) == pytest.approx(expected, rel=1e-6)
+
+
+def test_sidespin_does_not_deflect_a_rolling_ball() -> None:
+ """A rolling rigid sphere feels no lateral force from vertical-axis spin."""
+ table = TableParams(mu_roll=0.03, **SMOOTH_CLOTH)
+ ball = _launch(2.0, english_y=0.4, english_x=0.5)
+ assert abs(float(ball.omega[2])) > 1.0
+ stopped, _ = _roll_out(ball, table)
+ assert float(stopped.pos[1]) == pytest.approx(0.635, abs=1e-9)
+
+
+def test_ball_ball_collision_conserves_linear_momentum() -> None:
+ """Frictional collisions must conserve momentum to machine precision."""
+ params = BallParams()
+ table = TableParams()
+ a = Ball(
+ id=0,
+ pos=np.array([0.0, 0.0]),
+ vel=np.array([3.0, 0.7]),
+ omega=np.array([5.0, 60.0, -30.0]),
+ params=params,
+ )
+ b = Ball(
+ id=1,
+ pos=np.array([2 * params.radius, 0.0]),
+ vel=np.zeros(2),
+ omega=np.zeros(3),
+ params=params,
+ )
+ before = params.mass * (a.vel + b.vel)
+ a_out, b_out = resolve_ball_ball(a, b, table)
+ after = params.mass * (a_out.vel + b_out.vel)
+ np.testing.assert_allclose(after, before, atol=1e-12)
+
+
+def _kinetic_energy(balls: list[Ball]) -> float:
+ return float(
+ sum(
+ 0.5 * b.m * b.vel @ b.vel + 0.5 * b.I * b.omega @ b.omega
+ for b in balls
+ )
+ )
+
+
+def test_collisions_never_create_energy() -> None:
+ """Inelastic contacts are strictly dissipative, for ball-ball and cushions."""
+ params = BallParams()
+ table = TableParams()
+ a = Ball(
+ id=0,
+ pos=np.array([1.0, 0.6]),
+ vel=np.array([4.0, -1.2]),
+ omega=np.array([10.0, 80.0, 40.0]),
+ params=params,
+ )
+ b = Ball(
+ id=1,
+ pos=np.array([1.0 + 2 * params.radius, 0.6]),
+ vel=np.array([0.2, 0.0]),
+ omega=np.zeros(3),
+ params=params,
+ )
+ before = _kinetic_energy([a, b])
+ a_out, b_out = resolve_ball_ball(a, b, table)
+ assert _kinetic_energy([a_out, b_out]) <= before + 1e-12
+
+ rail = Ball(
+ id=2,
+ pos=np.array([params.radius - 1e-4, 0.6]),
+ vel=np.array([-3.0, 0.5]),
+ omega=np.array([0.0, 0.0, 25.0]),
+ params=params,
+ )
+ before_rail = _kinetic_energy([rail])
+ assert _kinetic_energy([resolve_cushion(rail, table)]) <= before_rail + 1e-12
+
+
+@pytest.mark.parametrize("speed", [1.0, 3.0])
+def test_cushion_restitution_matches_configured_coefficient(speed: float) -> None:
+ """Head-on rail rebound speed equals e_cushion x approach speed."""
+ params = BallParams()
+ table = TableParams(e_cushion=0.85)
+ ball = Ball(
+ id=0,
+ pos=np.array([params.radius - 1e-5, 0.635]),
+ vel=np.array([-speed, 0.0]),
+ omega=np.zeros(3),
+ params=params,
+ )
+ out = resolve_cushion(ball, table)
+ # Restitution is softened above 1.5 m/s to mimic cushion compliance.
+ expected = min(table.e_cushion, max(0.55, table.e_cushion - 0.03 * max(0.0, speed - 1.5)))
+ assert float(out.vel[0]) == pytest.approx(expected * speed, rel=1e-9)
+
+
+def test_draw_stun_follow_are_correctly_ordered() -> None:
+ """
+ Cue ball behaviour after a full-ball contact must follow the known ordering:
+ backspin draws it back behind contact, topspin sends it furthest forward.
+ """
+ sim = Simulator(dt=2e-4, max_time=12.0, collision_passes=10)
+ sim.table = TableParams(**SMOOTH_CLOTH)
+ contact_x = 1.2 - 2 * sim.ball_params.radius
+
+ endpoints = {}
+ for label, english_y in (("draw", -0.45), ("stun", 0.0), ("follow", 0.45)):
+ result = sim.simulate_shot(
+ ShotParams(speed=2.5, angle=0.0, english_y=english_y),
+ cue_pos=(0.6, 0.635),
+ obj_pos=(1.2, 0.635),
+ full_rack=False,
+ )
+ endpoints[label] = float(result.endpoints[0][0])
+
+ assert endpoints["draw"] < contact_x, "backspin must pull the cue ball back"
+ assert endpoints["draw"] < endpoints["stun"] < endpoints["follow"]
+
+
+@pytest.mark.parametrize(
+ "label,angle_deg,obj_pos,tolerance_m",
+ [
+ ("full ball", 0.0, (1.4, 0.635), 0.01),
+ # A thin cut amplifies the contact geometry, so it is the worst case for
+ # the timestep and gets the looser bound.
+ ("half ball", 2.0, (1.5, 0.66), 0.06),
+ ],
+)
+def test_training_timestep_is_converged(
+ label: str, angle_deg: float, obj_pos: tuple[float, float], tolerance_m: float
+) -> None:
+ """
+ Labels are generated at 2 ms, so halving the step must barely move the answer.
+
+ Without this bound the reported model errors could be measuring the
+ integrator's discretisation rather than the model, and no amount of training
+ would fix it. See docs/VALIDATION.md for the measured values.
+ """
+ table = TableParams(**SMOOTH_CLOTH)
+ shot = ShotParams(speed=3.0 if angle_deg else 2.5, angle=float(np.radians(angle_deg)))
+ cue_pos = (0.6, 0.60) if angle_deg else (0.6, 0.635)
+
+ endpoints = {}
+ for dt in (2e-3, 1e-3):
+ result = Simulator(table=table, dt=dt, max_time=8.0).simulate_shot(
+ shot, cue_pos=cue_pos, obj_pos=obj_pos
+ )
+ endpoints[dt] = (result.endpoints[0].copy(), result.endpoints[1].copy())
+ assert result.collision_events > 0, "the shot must actually make contact"
+
+ for index, ball in enumerate(("cue", "object")):
+ gap = float(np.linalg.norm(endpoints[2e-3][index] - endpoints[1e-3][index]))
+ assert gap < tolerance_m, f"{label} {ball} ball moved {gap * 1000:.0f} mm"
+
+
+def test_ghost_ball_geometry_pots_and_half_a_degree_misses() -> None:
+ """
+ Contact geometry against its closed form, and the tolerance that follows.
+
+ Aiming at the ghost-ball point — the object ball's centre pulled one diameter
+ back along the line to the pocket — must pot the ball. Half a degree off must
+ not, which is the measurement behind the claim that potting needs finer aim
+ than the surrogate can resolve.
+ """
+ table = TableParams(**SMOOTH_CLOTH)
+ sim = Simulator(table=table, dt=2e-3, max_time=8.0)
+ cue_pos, obj_pos = (0.60, 0.35), (1.70, 0.75)
+ pocket = np.array([table.length, table.width])
+
+ obj = np.asarray(obj_pos)
+ to_pocket = (pocket - obj) / np.linalg.norm(pocket - obj)
+ ghost = obj - 2 * sim.ball_params.radius * to_pocket
+ aim = float(np.arctan2(*(ghost - np.asarray(cue_pos))[::-1]))
+
+ def pots(angle: float) -> bool:
+ result = sim.simulate_shot(
+ ShotParams(speed=2.5, angle=angle), cue_pos=cue_pos, obj_pos=obj_pos
+ )
+ return bool(result.pocketed[1]) and not result.pocketed[0]
+
+ assert pots(aim), "the closed-form ghost-ball line must pot the object ball"
+ for offset_deg in (-0.5, 0.5):
+ assert not pots(aim + float(np.radians(offset_deg))), (
+ f"{offset_deg:+.1f}° off the ghost-ball line should miss"
+ )
+
+
+def test_analytic_baseline_tracks_simulator_on_open_table() -> None:
+ """The closed-form baseline agrees with the simulator when no rail is hit."""
+ table = TableParams(mu_slide=0.2, mu_roll=0.05, **SMOOTH_CLOTH)
+ sim = Simulator(table=table, dt=5e-4, max_time=20.0)
+ shot = ShotParams(speed=1.2, angle=0.0)
+ result = sim.simulate_shot(shot, cue_pos=(0.3, 0.635), object_ball=False)
+ predicted = analytic.predict_endpoint(shot, (0.3, 0.635), table)
+ assert float(result.endpoints[0][0]) == pytest.approx(float(predicted[0]), abs=0.02)
+
+
+def test_analytic_baseline_tracks_simulator_across_a_cushion() -> None:
+ """
+ Agreement must survive a rail, which is only true if the closed-form model
+ carries spin through the bounce: a ball arrives at the next rail sliding.
+ """
+ table = TableParams(mu_slide=0.2, mu_roll=0.02, **SMOOTH_CLOTH)
+ sim = Simulator(table=table, dt=5e-4, max_time=30.0)
+ shot = ShotParams(speed=2.0, angle=0.0)
+ start = (0.4, 0.5)
+ result = sim.simulate_shot(shot, cue_pos=start, object_ball=False)
+ assert result.cushion_events >= 1
+ predicted = analytic.predict_endpoint(shot, start, table)
+ error = float(np.linalg.norm(result.endpoints[0] - predicted))
+ assert error < 0.10, f"closed-form baseline off by {error * 1000:.0f} mm"
+
+
+def test_vertical_spin_decays_to_zero_without_chattering() -> None:
+ """
+ Spinning friction must land on zero, not step past it.
+
+ The decrement per timestep is constant, so subtracting it unconditionally
+ overshoots and flips the sign, leaving the ball spinning at a small
+ alternating rate for ever. Nothing downstream notices except that the table
+ never comes to rest, which is why this is asserted directly.
+ """
+ table = TableParams(**SMOOTH_CLOTH)
+ params = BallParams()
+ dt = 1e-3
+ # Deliberately start below one step's worth of decay, where overshoot bites.
+ step = 2.5 * table.mu_spin * G / params.radius * dt
+ ball = Ball(
+ id=0,
+ pos=np.array([1.0, 0.635]),
+ vel=np.zeros(2),
+ omega=np.array([0.0, 0.0, 0.4 * step]),
+ params=params,
+ )
+
+ previous = abs(float(ball.omega[2]))
+ for _ in range(50):
+ ball = integrate_ball(ball, table, dt)
+ current = abs(float(ball.omega[2]))
+ assert current <= previous + 1e-12, "spin grew while friction was removing it"
+ previous = current
+
+ assert float(ball.omega[2]) == 0.0
+ assert ball.motion_state(table) is MotionState.STATIONARY
+
+
+@pytest.mark.parametrize(
+ ("label", "shot"),
+ [
+ ("soft into the rack", ShotParams(speed=2.5, angle=0.12, english_x=0.2, english_y=-0.1)),
+ ("hard break", ShotParams(speed=8.0, angle=0.0)),
+ ("heavy right english", ShotParams(speed=4.0, angle=0.35, english_x=0.48)),
+ ("draw into the pack", ShotParams(speed=5.0, angle=-0.05, english_y=-0.45)),
+ ],
+)
+def test_every_shot_reaches_rest_before_the_time_limit(label: str, shot: ShotParams) -> None:
+ """
+ A shot must end because the balls stopped, not because the clock ran out.
+
+ Hitting ``max_time`` silently truncates the shot, and a truncated shot is
+ still reported as a resting position, so the failure is invisible in the
+ output and shows up only as a simulator that will not terminate.
+ """
+ sim = Simulator(dt=1e-3, max_time=15.0)
+ result = sim.simulate_shot(shot, full_rack=True, seed=7)
+ settled_at = float(result.times[-1])
+ assert settled_at < sim.max_time - 1.0, (
+ f"{label}: still moving at {settled_at:.1f}s, the {sim.max_time:.0f}s limit truncated it"
+ )
+
+
+def test_every_contact_in_a_rack_is_inside_the_contact_band() -> None:
+ """
+ A racked ball touches its neighbours, and the solver has to agree.
+
+ This is the regression test for a bug that was invisible in every other
+ check. The contact band used to be ``1e-4`` and the rack was built with a
+ ``1e-4`` clearance, so all thirty contacts in the triangle sat exactly on
+ the threshold that decides whether a contact exists. Which side of it each
+ one landed on came down to whether ``hypot`` rounded up or down: sixteen
+ registered, fourteen did not. A break then propagated through a contact
+ graph with holes in it, and balls in the middle of the rack came out of a
+ full-power break having never moved.
+
+ Neither the closed-form checks nor the parity harness could see it: the
+ single-ball mechanics were untouched, and the JavaScript port reproduced
+ the broken graph exactly, because it was a faithful port of it.
+ """
+ balls = make_full_rack(TableParams())[1:]
+ touching = 0
+ for i, a in enumerate(balls):
+ for b in balls[i + 1 :]:
+ gap = float(np.hypot(*(b.pos - a.pos))) - a.R - b.R
+ if gap > 1e-3:
+ continue # not neighbours in the triangle
+ touching += 1
+ assert gap <= CONTACT_BAND, (
+ f"balls {a.number} and {b.number} are {gap * 1e6:.3f} µm apart, "
+ f"outside the {CONTACT_BAND * 1e6:.0f} µm contact band, so the "
+ "solver will not see them as touching"
+ )
+ assert touching == 30, f"a five-row triangle has 30 contacts, found {touching}"
+
+
+def test_the_contact_band_is_far_from_both_things_that_could_swallow_it() -> None:
+ """
+ The band has to be wide against float noise and narrow against physics.
+
+ Stated as a test because the failure mode is silent either way: too tight
+ and real contacts are missed, too loose and balls collide with thin air.
+ """
+ noise = 1e-15 # the scale of rounding in a position, in metres
+ assert CONTACT_BAND / noise > 1000
+ assert BallParams().radius / CONTACT_BAND > 1000
+
+
+def test_resolving_a_resting_rack_changes_nothing_and_stops_immediately() -> None:
+ """
+ Thirty resting contacts must not cost thirty passes of work every step.
+
+ The sweep exits when a pass changes nothing, so a table that is merely
+ sitting there costs one pass. Without that, making the rack's contacts
+ visible would have quietly multiplied the cost of every timestep.
+ """
+ table = TableParams()
+ before = make_full_rack(table)
+ after = resolve_all_ball_collisions(before, table)
+ for a, b in zip(before, after):
+ assert np.allclose(a.pos, b.pos), f"ball {a.number} was moved while at rest"
+ assert np.allclose(a.vel, b.vel), f"ball {a.number} was pushed while at rest"
+
+
+def test_a_break_puts_every_ball_in_the_rack_in_motion() -> None:
+ """
+ The impulse has to reach the whole triangle, not the near half of it.
+
+ Two deliberately loose claims, because a break is chaotic and a tight
+ assertion about it is a lottery: every ball is set moving, and most of the
+ rack ends up somewhere else. What the numbers are is not the point; that
+ there is no ball the impulse never reached is.
+
+ The distribution across those balls is much less even than a real break —
+ see the multi-contact note in ``docs/VALIDATION.md``. Resolving contacts
+ pairwise in sequence is an approximation to a rack in which fifteen balls
+ are touching at once, and it is the largest known departure from reality in
+ this simulator.
+ """
+ sim = Simulator(dt=1e-3, max_time=20.0)
+ result = sim.simulate_shot(ShotParams(speed=8.0, angle=0.0), full_rack=True, seed=7)
+
+ untouched, displaced = [], 0
+ for ball_id, path in result.trajectories.items():
+ if ball_id == 0:
+ continue
+ if float(np.max(np.linalg.norm(result.velocities[ball_id], axis=1))) < 1e-3:
+ untouched.append(ball_id)
+ displaced += float(np.linalg.norm(path[-1] - path[0])) > 0.05
+
+ assert not untouched, f"balls {sorted(untouched)} never moved on a full-power break"
+ assert displaced >= 8, f"only {displaced} of 15 balls left the rack area"
+
+
+def test_a_harder_break_opens_the_table_further() -> None:
+ """
+ The property whose sign was wrong, asserted so it cannot go wrong quietly.
+
+ Nobody writes this test first, because it is too obvious to be worth
+ stating. That is exactly why the contact-band defect survived: every
+ closed-form check passed, the browser port agreed to eleven decimal places,
+ and the only visible symptom was that a rack struck at 8 m/s ended up
+ *tighter* than one struck at 5. The measurement that found it is the one
+ below.
+
+ Coarse on purpose. A 3 ms step and one rack per speed keeps this to a few
+ seconds, and the claim is about the direction of an aggregate, which is not
+ a quantity a millimetre of discretisation decides.
+ """
+ spreads = []
+ for speed in (3.0, 6.0, 9.0):
+ sim = Simulator(dt=3e-3, max_time=15.0)
+ result = sim.simulate_shot(ShotParams(speed=speed, angle=0.0), full_rack=True, seed=7)
+ resting = np.array(
+ [
+ result.endpoints[ball_id][:2]
+ for ball_id in result.trajectories
+ if ball_id != 0 and not result.pocketed[ball_id]
+ ]
+ )
+ centre = resting.mean(axis=0)
+ spreads.append(float(np.mean(np.linalg.norm(resting - centre, axis=1))))
+
+ assert spreads[0] < spreads[1] < spreads[2], (
+ "a harder break has to leave the balls further apart, but the mean "
+ f"distance from the centre of the pack went {spreads[0]:.3f} → "
+ f"{spreads[1]:.3f} → {spreads[2]:.3f} m at 3, 6 and 9 m/s"
+ )
diff --git a/web/css/style.css b/web/css/style.css
new file mode 100644
index 0000000..2e98a4e
--- /dev/null
+++ b/web/css/style.css
@@ -0,0 +1,735 @@
+:root {
+ --bg: #0b0f14;
+ --panel: #131a22;
+ --panel-2: #1a232d;
+ --line: #24303c;
+ --ink: #e6edf3;
+ --ink-dim: #8b9bab;
+ --ink-faint: #5d6b7a;
+ --accent: #4aa8ff;
+ --good: #3fb950;
+ --warn: #f0883e;
+ --bad: #f85149;
+ --mono: ui-monospace, "SF Mono", "Cascadia Mono", Menlo, Consolas, monospace;
+ --sans: system-ui, -apple-system, "Segoe UI", Inter, Roboto, sans-serif;
+}
+
+* {
+ box-sizing: border-box;
+}
+
+html,
+body {
+ margin: 0;
+ padding: 0;
+ background: var(--bg);
+ color: var(--ink);
+ font-family: var(--sans);
+ -webkit-font-smoothing: antialiased;
+}
+
+a {
+ color: var(--accent);
+ text-decoration: none;
+}
+a:hover {
+ text-decoration: underline;
+}
+
+.wrap {
+ max-width: 1400px;
+ margin: 0 auto;
+ padding: 22px 20px 80px;
+}
+
+/* ---------- header ---------- */
+
+header.masthead {
+ display: flex;
+ align-items: baseline;
+ gap: 16px;
+ flex-wrap: wrap;
+ border-bottom: 1px solid var(--line);
+ padding-bottom: 14px;
+ margin-bottom: 20px;
+}
+
+.masthead h1 {
+ font-size: 21px;
+ margin: 0;
+ letter-spacing: -0.2px;
+}
+
+.masthead h1 span {
+ color: var(--accent);
+}
+
+.masthead .tag {
+ color: var(--ink-dim);
+ font-size: 14px;
+ flex: 1;
+ min-width: 260px;
+}
+
+.masthead nav {
+ display: flex;
+ gap: 14px;
+ font-size: 13px;
+}
+
+/* ---------- layout ---------- */
+
+.stage {
+ display: grid;
+ grid-template-columns: minmax(0, 1fr) 330px;
+ gap: 20px;
+ align-items: start;
+}
+
+@media (max-width: 1080px) {
+ .stage {
+ grid-template-columns: minmax(0, 1fr);
+ }
+}
+
+.table-column {
+ min-width: 0;
+}
+
+/* ---------- table ---------- */
+
+.table-shell {
+ position: relative;
+ border-radius: 14px;
+ overflow: hidden;
+ background: #0a0d11;
+ box-shadow: 0 20px 50px rgba(0, 0, 0, 0.55);
+}
+
+#table {
+ display: block;
+ width: 100%;
+ height: auto;
+ cursor: crosshair;
+ touch-action: none;
+}
+
+.banner {
+ position: absolute;
+ inset: 0;
+ display: none;
+ align-items: center;
+ justify-content: center;
+ flex-direction: column;
+ gap: 14px;
+ background: rgba(6, 10, 14, 0.86);
+ backdrop-filter: blur(3px);
+ text-align: center;
+ padding: 20px;
+}
+
+.banner.show {
+ display: flex;
+}
+
+.banner h2 {
+ margin: 0;
+ font-size: 30px;
+ letter-spacing: -0.5px;
+}
+
+.banner p {
+ margin: 0;
+ color: var(--ink-dim);
+ max-width: 420px;
+ line-height: 1.55;
+}
+
+/* ---------- controls ---------- */
+
+.controls {
+ display: grid;
+ grid-template-columns: auto 1fr auto;
+ gap: 18px;
+ align-items: center;
+ background: var(--panel);
+ border: 1px solid var(--line);
+ border-top: none;
+ border-radius: 0 0 14px 14px;
+ padding: 14px 16px;
+ margin-top: -14px;
+ padding-top: 24px;
+}
+
+@media (max-width: 620px) {
+ .controls {
+ grid-template-columns: 1fr;
+ }
+}
+
+/* ---------- shot log ---------- */
+
+.shotlog {
+ margin-top: 14px;
+ background: var(--panel);
+ border: 1px solid var(--line);
+ border-radius: 12px;
+ padding: 12px 14px 6px;
+}
+
+.shotlog-head {
+ display: flex;
+ justify-content: space-between;
+ align-items: baseline;
+ gap: 12px;
+ font-size: 10.5px;
+ text-transform: uppercase;
+ letter-spacing: 0.7px;
+ color: var(--ink-dim);
+}
+
+.shotlog-head span + span {
+ text-transform: none;
+ letter-spacing: 0;
+ font-size: 11.5px;
+ color: var(--ink-faint);
+}
+
+.shotlog ol {
+ list-style: none;
+ margin: 8px 0 0;
+ padding: 0;
+ font-family: var(--mono);
+ font-size: 11.5px;
+ line-height: 1.55;
+}
+
+.shotlog li {
+ display: flex;
+ gap: 10px;
+ padding: 5px 0;
+ border-top: 1px solid var(--line);
+ color: var(--ink-dim);
+}
+
+.shotlog li:first-child {
+ border-top: none;
+}
+
+.shotlog li.empty {
+ color: var(--ink-faint);
+ font-family: inherit;
+}
+
+.shotlog .who {
+ flex: 0 0 34px;
+ font-weight: 700;
+ color: var(--ink);
+}
+
+.shotlog .who.bot {
+ color: var(--warn);
+}
+
+.shotlog .what {
+ flex: 1 1 auto;
+ min-width: 0;
+}
+
+.shotlog .why {
+ flex: 0 0 auto;
+ color: var(--ink-faint);
+}
+
+.shotlog .why.foul {
+ color: var(--warn);
+}
+
+.spin-widget {
+ text-align: center;
+}
+
+.spin-widget canvas {
+ display: block;
+ cursor: pointer;
+ border-radius: 50%;
+}
+
+.power-block label,
+.spin-widget .caption {
+ display: block;
+ font-size: 10.5px;
+ text-transform: uppercase;
+ letter-spacing: 0.9px;
+ color: var(--ink-faint);
+ margin-bottom: 7px;
+}
+
+.power-row {
+ display: flex;
+ align-items: center;
+ gap: 12px;
+}
+
+input[type="range"] {
+ -webkit-appearance: none;
+ appearance: none;
+ flex: 1;
+ height: 6px;
+ border-radius: 3px;
+ background: linear-gradient(90deg, #2b6cb0, #f0883e, #f85149);
+ outline: none;
+}
+
+input[type="range"]::-webkit-slider-thumb {
+ -webkit-appearance: none;
+ width: 18px;
+ height: 18px;
+ border-radius: 50%;
+ background: var(--ink);
+ border: 2px solid var(--bg);
+ cursor: pointer;
+}
+
+input[type="range"]::-moz-range-thumb {
+ width: 16px;
+ height: 16px;
+ border-radius: 50%;
+ background: var(--ink);
+ border: 2px solid var(--bg);
+ cursor: pointer;
+}
+
+.readout {
+ font-family: var(--mono);
+ font-size: 13px;
+ color: var(--ink);
+ min-width: 74px;
+ text-align: right;
+}
+
+button {
+ font-family: var(--sans);
+ font-size: 13px;
+ font-weight: 600;
+ color: var(--ink);
+ background: var(--panel-2);
+ border: 1px solid var(--line);
+ border-radius: 8px;
+ padding: 9px 15px;
+ cursor: pointer;
+ transition: border-color 0.12s, background 0.12s;
+}
+
+button:hover:not(:disabled) {
+ border-color: var(--accent);
+ background: #202b37;
+}
+
+button:disabled {
+ opacity: 0.4;
+ cursor: not-allowed;
+}
+
+button.primary {
+ background: #1b4f80;
+ border-color: #2d6ba8;
+}
+
+button.primary:hover:not(:disabled) {
+ background: #226199;
+}
+
+.button-stack {
+ display: flex;
+ flex-direction: column;
+ gap: 8px;
+}
+
+/* ---------- side panels ---------- */
+
+.panel {
+ background: var(--panel);
+ border: 1px solid var(--line);
+ border-radius: 12px;
+ padding: 14px 15px;
+ margin-bottom: 14px;
+}
+
+.panel h3 {
+ margin: 0 0 11px;
+ font-size: 10.5px;
+ text-transform: uppercase;
+ letter-spacing: 1px;
+ color: var(--ink-faint);
+ font-weight: 700;
+}
+
+.turn-line {
+ display: flex;
+ align-items: center;
+ gap: 9px;
+ font-size: 16px;
+ font-weight: 650;
+ margin-bottom: 6px;
+}
+
+.dot {
+ width: 9px;
+ height: 9px;
+ border-radius: 50%;
+ background: var(--accent);
+ flex: none;
+}
+
+.dot.bot {
+ background: var(--warn);
+}
+
+.message {
+ font-size: 13px;
+ color: var(--ink-dim);
+ line-height: 1.5;
+ min-height: 38px;
+}
+
+.message .foul {
+ color: var(--bad);
+ font-weight: 600;
+}
+
+.groups {
+ display: grid;
+ grid-template-columns: 1fr 1fr;
+ gap: 10px;
+ margin-top: 12px;
+}
+
+.group-card {
+ background: var(--panel-2);
+ border: 1px solid var(--line);
+ border-radius: 8px;
+ padding: 9px 10px;
+}
+
+.group-card.active {
+ border-color: var(--accent);
+}
+
+.group-card .who {
+ font-size: 10.5px;
+ text-transform: uppercase;
+ letter-spacing: 0.8px;
+ color: var(--ink-faint);
+}
+
+.group-card .what {
+ font-size: 14px;
+ font-weight: 650;
+ margin-top: 3px;
+}
+
+.ball-row {
+ display: flex;
+ flex-wrap: wrap;
+ gap: 4px;
+ margin-top: 7px;
+ min-height: 18px;
+}
+
+.pip {
+ width: 18px;
+ height: 18px;
+ border-radius: 50%;
+ font-size: 9.5px;
+ font-weight: 700;
+ color: #101010;
+ display: flex;
+ align-items: center;
+ justify-content: center;
+ border: 1px solid rgba(0, 0, 0, 0.45);
+}
+
+.pip.stripe {
+ background-image: linear-gradient(
+ to bottom,
+ #f2f2ec 0 26%,
+ var(--pip) 26% 74%,
+ #f2f2ec 74% 100%
+ );
+}
+
+/* ---------- inspector ---------- */
+
+.stat-grid {
+ display: grid;
+ grid-template-columns: 1fr 1fr;
+ gap: 8px 12px;
+ font-family: var(--mono);
+ font-size: 12px;
+}
+
+.stat {
+ display: flex;
+ justify-content: space-between;
+ gap: 8px;
+ border-bottom: 1px dotted var(--line);
+ padding-bottom: 3px;
+}
+
+.stat .k {
+ color: var(--ink-faint);
+}
+
+.stat .v {
+ color: var(--ink);
+ font-variant-numeric: tabular-nums;
+}
+
+.state-chip {
+ display: inline-block;
+ font-family: var(--mono);
+ font-size: 11px;
+ font-weight: 700;
+ letter-spacing: 0.5px;
+ padding: 3px 9px;
+ border-radius: 5px;
+ background: var(--panel-2);
+ border: 1px solid var(--line);
+ color: var(--ink-dim);
+ margin-bottom: 10px;
+}
+
+.state-chip.sliding {
+ color: #ffd479;
+ border-color: #6b5322;
+ background: #2a2213;
+}
+.state-chip.rolling {
+ color: var(--good);
+ border-color: #26562f;
+ background: #13251a;
+}
+.state-chip.spinning {
+ color: #d2a8ff;
+ border-color: #4a3568;
+ background: #221a2e;
+}
+.state-chip.stationary {
+ color: var(--ink-faint);
+}
+.state-chip.pocketed {
+ color: var(--warn);
+ border-color: #6b3a22;
+ background: #2a1a13;
+}
+
+#trace {
+ display: block;
+ width: 100%;
+ height: auto;
+ margin-top: 11px;
+ border-radius: 6px;
+ background: #0d1319;
+}
+
+.trace-legend {
+ display: flex;
+ gap: 13px;
+ font-size: 10.5px;
+ color: var(--ink-faint);
+ margin-top: 6px;
+ flex-wrap: wrap;
+}
+
+.trace-legend i {
+ display: inline-block;
+ width: 9px;
+ height: 2px;
+ vertical-align: middle;
+ margin-right: 4px;
+}
+
+/* ---------- bot panel ---------- */
+
+.bot-line {
+ font-family: var(--mono);
+ font-size: 11.5px;
+ color: var(--ink-dim);
+ line-height: 1.7;
+}
+
+.bot-line b {
+ color: var(--ink);
+ font-weight: 600;
+}
+
+.progress {
+ height: 3px;
+ border-radius: 2px;
+ background: var(--panel-2);
+ overflow: hidden;
+ margin: 9px 0;
+}
+
+.progress > div {
+ height: 100%;
+ width: 0;
+ background: var(--warn);
+ transition: width 0.08s linear;
+}
+
+.setting-row {
+ display: flex;
+ align-items: center;
+ justify-content: space-between;
+ gap: 10px;
+ font-size: 12.5px;
+ color: var(--ink-dim);
+ padding: 5px 0;
+}
+
+select {
+ font-family: var(--sans);
+ font-size: 12.5px;
+ color: var(--ink);
+ background: var(--panel-2);
+ border: 1px solid var(--line);
+ border-radius: 6px;
+ padding: 5px 8px;
+ cursor: pointer;
+}
+
+input[type="checkbox"] {
+ accent-color: var(--accent);
+ cursor: pointer;
+ width: 15px;
+ height: 15px;
+}
+
+/* ---------- explainer ---------- */
+
+.explainer {
+ margin-top: 42px;
+ border-top: 1px solid var(--line);
+ padding-top: 30px;
+}
+
+.explainer > h2 {
+ font-size: 22px;
+ margin: 0 0 8px;
+ letter-spacing: -0.3px;
+}
+
+.explainer > p.lede {
+ color: var(--ink-dim);
+ font-size: 14.5px;
+ line-height: 1.65;
+ max-width: 760px;
+ margin: 0 0 26px;
+}
+
+.cards {
+ display: grid;
+ grid-template-columns: repeat(auto-fit, minmax(310px, 1fr));
+ gap: 16px;
+}
+
+.card {
+ background: var(--panel);
+ border: 1px solid var(--line);
+ border-radius: 12px;
+ padding: 17px 18px;
+}
+
+.card h3 {
+ margin: 0 0 9px;
+ font-size: 14.5px;
+ color: var(--ink);
+}
+
+.card p {
+ margin: 0 0 10px;
+ font-size: 13.2px;
+ line-height: 1.62;
+ color: var(--ink-dim);
+}
+
+.card p:last-child {
+ margin-bottom: 0;
+}
+
+.card code,
+.explainer code {
+ font-family: var(--mono);
+ font-size: 12px;
+ background: var(--panel-2);
+ border: 1px solid var(--line);
+ border-radius: 4px;
+ padding: 1px 5px;
+ color: #a5d6ff;
+}
+
+.eq {
+ display: block;
+ font-family: var(--mono);
+ font-size: 12.5px;
+ color: #a5d6ff;
+ background: #0d1319;
+ border: 1px solid var(--line);
+ border-left: 2px solid var(--accent);
+ border-radius: 5px;
+ padding: 9px 12px;
+ margin: 11px 0;
+ overflow-x: auto;
+}
+
+.measure {
+ display: flex;
+ flex-wrap: wrap;
+ gap: 10px;
+ margin-top: 12px;
+}
+
+.measure div {
+ background: var(--panel-2);
+ border: 1px solid var(--line);
+ border-radius: 8px;
+ padding: 8px 12px;
+ flex: 1 1 120px;
+}
+
+.measure .n {
+ font-family: var(--mono);
+ font-size: 16px;
+ font-weight: 600;
+ color: var(--good);
+}
+
+.measure .l {
+ font-size: 10.5px;
+ color: var(--ink-faint);
+ text-transform: uppercase;
+ letter-spacing: 0.6px;
+ margin-top: 2px;
+}
+
+footer {
+ margin-top: 42px;
+ padding-top: 18px;
+ border-top: 1px solid var(--line);
+ font-size: 12.5px;
+ color: var(--ink-faint);
+ line-height: 1.7;
+}
+
+.kbd {
+ font-family: var(--mono);
+ font-size: 11px;
+ background: var(--panel-2);
+ border: 1px solid var(--line);
+ border-bottom-width: 2px;
+ border-radius: 4px;
+ padding: 1px 5px;
+ color: var(--ink-dim);
+}
diff --git a/web/data/facts.json b/web/data/facts.json
new file mode 100644
index 0000000..1601b90
--- /dev/null
+++ b/web/data/facts.json
@@ -0,0 +1,18 @@
+{
+ "physics-tests": "32",
+ "parity-worst": "0.0012 mm",
+ "parity-cases": "35",
+ "parity-speedup": "65\u00d7",
+ "parity-table-seconds": "156 s",
+ "parity-python-seconds": "39 s",
+ "parity-browser-seconds": "0.59 s",
+ "surrogate-error": "376 mm",
+ "analytic-error": "494 mm",
+ "gbm-error": "382 mm",
+ "samples": "20,000",
+ "contact-cuenet": "696 mm",
+ "contact-gbm": "610 mm",
+ "sim-full-rack": "4.6 s",
+ "surrogate-latency": "0.60 ms",
+ "closed-form-latency": "0.27 ms"
+}
diff --git a/web/data/parity.json b/web/data/parity.json
new file mode 100644
index 0000000..b832b5b
--- /dev/null
+++ b/web/data/parity.json
@@ -0,0 +1,10 @@
+{
+ "cases": 35,
+ "agreed": 35,
+ "worst_mm": 0.0011522116664723693,
+ "worst_case": "break-6ms",
+ "table_seconds": 155.76400000000493,
+ "python_seconds": 38.60540390199458,
+ "browser_seconds": 0.5890029370000001,
+ "speedup": 65.54365263206586
+}
diff --git a/web/index.html b/web/index.html
new file mode 100644
index 0000000..43389f6
--- /dev/null
+++ b/web/index.html
@@ -0,0 +1,349 @@
+
+
+
+ This is not a pool game with physics bolted on. It is a rigid-body simulator that + was written first, checked against closed-form mechanics, and then given a table + and an opponent so the model can be watched while it runs. Everything below is + measured by something in the repository, and the command that measures it is named. +
+ ++ A struck ball does not simply decelerate. While its contact point is slipping + against the cloth, friction acts opposite the slip velocity rather + than the ball's velocity, and the same force exerts a torque that spins the + ball up: +
+u = v + ω × (−R ẑ), a = −μs g û, α = (−R ẑ × F)/I
+
+ Slip decays 3.5 times faster than the centre of mass slows, so the ball settles
+ into rolling at exactly 5/7 of its launch speed and only then
+ switches to the much smaller rolling resistance. Watch the
+ Cue ball, live panel during a shot: the yellow slip curve collapses to
+ zero and the green speed curve lands on the blue line. That line is a
+ prediction from the mechanics, not a fitted parameter.
+
+ Simulators are easy to make self-consistent and wrong. The test suite compares
+ the integrator against results derived independently on paper, not against its
+ own earlier output: the 5/7 rolling speed, the
+ 12v₀²/(49 μs g) sliding distance, spin decay, cushion
+ restitution, the draw/stun/follow ordering, exact momentum conservation and
+ strictly decreasing energy.
+
+ An earlier version of this repository failed all of that. The contact-point
+ slip and the friction torque had opposite handedness, so friction drove the
+ ball away from rolling: measured rolling speed was zero instead of
+ 5/7 v₀, and stopping distance was four times short.
+
+ Closed-form checks exercise one ball at a time, so a defect that only exists
+ between balls survives them. This one did. The solver counted two balls as
+ touching when the gap between their surfaces was under
+ 1e-4 m — and the rack was built with a 1e-4 m
+ clearance. All thirty contacts in the triangle sat exactly on the threshold,
+ and which side each fell on came down to whether hypot rounded up
+ or down. Sixteen registered. Fourteen did not.
+
+ A break propagated through a contact graph with holes in it: balls in the + middle of the rack came out of a full-power break having barely moved, and the + table opened up less the harder it was struck. It was found by + measuring break spread against cue speed and getting the sign wrong, not by a + test failing. What fixed it was giving the tolerance a name, using the same one + everywhere, and racking the balls actually touching so nothing sits on the + boundary. A test now asserts the property that caught it, because it is + exactly the kind nobody writes down: a harder break opens the table further. +
++ Measuring that turned up something the fix does not solve, and which is worth + stating rather than hiding. A rack is resolved as a chain of about fifteen + pairwise collisions, and each one applies restitution, so the survival factor + compounds: only 48% of the kinetic energy comes through a 10 m/s break, + where the real event is one stress wave that dissipates once. The break you + see below opens the table, but less than a real one would. It is the largest + known departure from reality here, and fixing it means treating a simultaneous + contact set as a single event rather than adjusting a coefficient. +
+
+ web/js/physics.js is a hand port of
+ src/cueai/physics/, and a port nobody measures is a rumour.
+ scripts/export_parity_cases.py runs
+ 35 shots through the Python simulator —
+ draw, follow, english off two rails, thin cuts, clusters, and full sixteen-ball
+ breaks — and records where every ball stopped.
+ web/test/parity.mjs replays them here and compares.
+
+ The largest disagreement across all of them is + 0.0012 mm — a thousandth of a + millimetre, accumulated over seconds of a chaotic sixteen-ball break, and far + below the scale at which the two could be said to behave differently. It runs + in continuous integration, so they cannot quietly drift apart. +
++ Twenty headless bot-against-bot games walk every branch of the rules, and + check on each of the roughly eight hundred thousand physics steps they take + that no two balls are ever sharing space and nothing is ever inside a cushion. + The worst overlap seen is 0.5 mm on a 57 mm ball, a fifth of a pixel at + the size this table is drawn. They also assert that every break moves at + least ten of the fifteen — it moves twelve — because a break that clips the + apex leaves the rack standing and looks exactly like a limitation of the + physics rather than a badly aimed cue. +
++ Aiming a pot needs no learned model: the ghost-ball construction is exact, and + a test asserts that aiming at it drops the ball while half a degree either + side misses. So the bot spends nothing on aiming. It enumerates every + ball-and-pocket pair in closed form, discards whatever is blocked or cut too + thin, and spends its whole budget simulating the survivors to see + what each one leaves behind. +
++ That is the argument this repository makes about fast physics, running in + front of you: the value of a cheap simulation is the number of futures you can + afford to look at. The panel above reports how many candidates the bot + enumerated, how many it actually simulated, and how much table time that was. +
++ The Python side also trains a surrogate that predicts resting positions + directly, at 0.60 ms against + 4.6 s for a full rack in the reference + simulator. It cuts mean endpoint error from + 494 mm to + 376 mm against the closed-form + baseline it corrects, over 20,000 simulated + shots. +
++ The bot above does not use it, and the reason is the more useful result. Ported + to the browser the same physics runs + 65× faster, which is enough to simulate + the shots worth considering exactly. A surrogate is worth its error when you + need to screen far more candidates than you can afford to simulate; inside one + turn, you do not. Reported against interest: once a ball-ball collision has to + be modelled, it loses to plain gradient boosting + (696 mm against + 610 mm). Every number here is rewritten by + the training run. +
+Nothing here is a screenshot of a number that was true once.
+make check # ruff, mypy, the whole test suite
+make parity # export shots from Python, replay them in Node
+make selfplay # 20 headless games, every rule branch
+make browser # load this page in Chrome and play a game
+make all # dataset → training → benchmarks → figures
+
+ The figures on this page are written by
+ scripts/site_facts.py from those runs rather than typed in, and
+ continuous integration fails if they drift. The game itself is dependency-free
+ ES modules: no bundler, no framework, no build step.
+ python3 -m http.server in web/ is enough to run
+ everything you see.
+