Skip to content

Replace in-house quartic solver with 1010 algorithm#235

Merged
ekiefl merged 6 commits into
mainfrom
ek/1010-quartics
Nov 12, 2025
Merged

Replace in-house quartic solver with 1010 algorithm#235
ekiefl merged 6 commits into
mainfrom
ek/1010-quartics

Conversation

@ekiefl
Copy link
Copy Markdown
Owner

@ekiefl ekiefl commented Nov 9, 2025

Improvement in accuracy

error_cdf

Figure: Relative error of roots (abs(roots) - abs(true_roots)) / abs(true_roots) for 100,000 "difficult" coefficients.

Legend:

PR = polynomial-roots

OQS = 1010-algoritm

pooltool = the old solver for pooltool

Numba = the new solver for pooltool

Improvement in speed

Old: 336.566 ms / 100000 solves
New: 241.622 ms / 100000 solves

Additional speed tests: ~/Software/pooltool_ml/pooltool/tests/ptmath/roots/data/README.md:

Batch size: 10 equations
C (ctypes):    0.1270 ms ± 0.0551 ms
Numba:         0.0125 ms ± 0.0001 ms

Batch size: 100 equations
C (ctypes):    1.1824 ms ± 0.2586 ms
Numba:         0.1272 ms ± 0.0234 ms

Batch size: 1000 equations
C (ctypes):   12.0158 ms ± 8.6863 ms
Numba:         1.3218 ms ± 0.1459 ms

Batch size: 10000 equations
C (ctypes):  109.2857 ms ± 5.3903 ms
Numba:        11.7594 ms ± 0.4457 ms

Single roots: 10000 polynomials
C (ctypes):   11.0477 μs ± 15.3976 μs
Numba:         2.9295 μs ± 1.7345 μs
C (direct):    0.3394 μs ± 0.5790 μs

Main takeaways here:

  • ctypes is slow because of overhead
  • direct C implementation is crazy fast

Summary by CodeRabbit

  • Refactor

    • Removed configurable quartic_solver parameter from simulation and event detection APIs; solver behavior is now internally unified.
    • Optimized polynomial root-finding with Numba acceleration.
  • Tests

    • Introduced data-driven test suite with precomputed ground truth data for polynomial root verification.

@ekiefl
Copy link
Copy Markdown
Owner Author

ekiefl commented Nov 9, 2025

@CodeRabbit review

@coderabbitai
Copy link
Copy Markdown

coderabbitai Bot commented Nov 9, 2025

Walkthrough

This PR removes a pluggable QuarticSolver enum-based architecture throughout the codebase and replaces it with a single hard-coded Numba-based quartic solver. Function signatures accepting a solver parameter are simplified across simulation and root-solving modules, while a new NumPy/Numba solver implementation is introduced. Test infrastructure is updated to remove solver parameterization and employ data-driven validation.

Changes

Cohort / File(s) Change Summary
API simplification and solver removal
pooltool/ptmath/roots/quartic.py, pooltool/evolution/event_based/simulate.py, pooltool/evolution/event_based/introspection.py
Removed QuarticSolver enum class; eliminated solver parameter from solve_quartics(), simulate(), get_next_event(), collision detection functions, and simulate_with_snapshots(); simplified all call sites to remove solver arguments
New Numba-based solver implementation
pooltool/ptmath/roots/_quartic_numba.py
Introduced complete quartic root-solver with solve(), solve_many(), and helper functions using NumPy and Numba; includes Newton-Raphson refinement, edge-case handling, and numerical stability checks for batch processing
Test suite refactoring
tests/evolution/event_based/test_simulate.py, tests/ptmath/roots/test_quartic.py
Removed pytest parameterization and solver arguments from test functions; replaced parameterized tests with data-driven tests using solve_many() and root-matching validation against ground-truth data
Test data infrastructure
tests/ptmath/roots/_generate_quartic_coeffs.py, tests/ptmath/roots/README.md
Added data generation script with functions to produce quartic coefficients from specified root configurations using mpmath for high precision; documented ground-truth test data generation and storage format

Estimated code review effort

🎯 4 (Complex) | ⏱️ ~60 minutes

  • Critical file: pooltool/ptmath/roots/_quartic_numba.py — Dense numerical logic with multiple root-finding paths, Newton-Raphson refinement, fallback strategies, and edge-case handling; requires verification of mathematical correctness and numerical stability
  • Interdependency risk: Changes span core simulation pathways; must verify that all solver invocations are consistently updated and produce equivalent or superior numerical results
  • Test validation: New data-driven tests rely on high-precision ground-truth generation; verify test data generation logic and tolerance thresholds are appropriate for the new solver
  • Behavioral impact: Switching to a single solver implementation may affect collision detection accuracy or simulation performance; review test assertions and expected values for any divergence

Possibly related PRs

Pre-merge checks and finishing touches

✅ Passed checks (3 passed)
Check name Status Explanation
Title check ✅ Passed The title 'Replace in-house quartic solver with 1010 algorithm' accurately describes the main change: replacing the previous quartic solving pipeline with a new NumPy/Numba-based implementation (the 1010 algorithm), which is evident across all modified files.
Docstring Coverage ✅ Passed No functions found in the changed files to evaluate docstring coverage. Skipping docstring coverage check.
Description Check ✅ Passed Check skipped - CodeRabbit’s high-level summary is enabled.
✨ Finishing touches
  • 📝 Generate docstrings
🧪 Generate unit tests (beta)
  • Create PR with unit tests
  • Post copyable unit tests in a comment
  • Commit unit tests in branch ek/1010-quartics

Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out.

❤️ Share

Comment @coderabbitai help to get the list of available commands and usage tips.

@ekiefl
Copy link
Copy Markdown
Owner Author

ekiefl commented Nov 9, 2025

@coderabbitai review

@coderabbitai
Copy link
Copy Markdown

coderabbitai Bot commented Nov 9, 2025

✅ Actions performed

Review triggered.

Note: CodeRabbit is an incremental review system and does not re-review already reviewed commits. This command is applicable only when automatic reviews are paused.

1 similar comment
@coderabbitai
Copy link
Copy Markdown

coderabbitai Bot commented Nov 9, 2025

✅ Actions performed

Review triggered.

Note: CodeRabbit is an incremental review system and does not re-review already reviewed commits. This command is applicable only when automatic reviews are paused.

@codecov
Copy link
Copy Markdown

codecov Bot commented Nov 9, 2025

Codecov Report

❌ Patch coverage is 10.70111% with 484 lines in your changes missing coverage. Please review.
✅ Project coverage is 44.84%. Comparing base (0e89b7f) to head (7ab6d52).
⚠️ Report is 42 commits behind head on main.

Files with missing lines Patch % Lines
pooltool/ptmath/roots/_quartic_numba.py 6.91% 458 Missing ⚠️
pooltool/ptmath/roots/core.py 17.85% 23 Missing ⚠️
pooltool/evolution/event_based/solve.py 70.00% 3 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main     #235      +/-   ##
==========================================
+ Coverage   44.22%   44.84%   +0.61%     
==========================================
  Files         127      141      +14     
  Lines        8554     9767    +1213     
==========================================
+ Hits         3783     4380     +597     
- Misses       4771     5387     +616     
Flag Coverage Δ
service 44.84% <10.70%> (+0.61%) ⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

Copy link
Copy Markdown

@coderabbitai coderabbitai Bot left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actionable comments posted: 2

Caution

Some comments are outside the diff and can’t be posted inline due to platform limitations.

⚠️ Outside diff range comments (1)
pooltool/ptmath/roots/quartic.py (1)

12-30: Handle degenerate polynomials before calling the Numba solver

The new _solver returns [0, 0, 0, 0] when the leading quartic term is zero, which propagates through get_real_positive_smallest_roots as “no positive root”. Any legitimate solution of the remaining cubic/quadratic therefore gets discarded, and downstream collision code skips real events. The old pluggable solvers fell back to lower-degree formulas, so this is a regression.

We can keep the fast path for true quartics while routing degenerate rows through np.roots (or similar) before taking real-positive minima:

 def solve_quartics(ps: NDArray[np.float64]) -> NDArray[np.float64]:
-    """Returns the smallest positive and real root for each quartic polynomial.
-
-    Args:
-        ps:
-            A mx5 array of polynomial coefficients, where m is the number of equations.
-            The columns are in the order a, b, c, d, e, where these coefficients make up
-            the quartic polynomial equation at^4 + bt^3 + ct^2 + dt + e = 0.
-        solver:
-            The method used to calculate the roots. See
-            pooltool.ptmath.roots.quartic.QuarticSolver.
-
-    Returns:
-        roots:
-            An array of shape m. Each value is the smallest root that is real and
-            positive. If no such root exists (e.g. all roots have complex), then
-            `np.inf` is returned.
-    """
-    return get_real_positive_smallest_roots(_solver(ps))
+    """Returns the smallest positive and real root for each polynomial."""
+    ps = np.asarray(ps, dtype=np.float64)
+    ps = np.atleast_2d(ps)
+
+    leading_zero = np.isclose(ps[:, 0], 0.0)
+    if not leading_zero.any():
+        return get_real_positive_smallest_roots(_solver(ps))
+
+    all_roots = np.full((ps.shape[0], 4), np.nan + 0j, dtype=np.complex128)
+
+    if (~leading_zero).any():
+        all_roots[~leading_zero] = _solver(ps[~leading_zero])
+
+    for idx in np.nonzero(leading_zero)[0]:
+        trimmed = np.trim_zeros(ps[idx], trim="f")
+        if trimmed.size <= 1:
+            continue  # constant polynomial -> no finite roots
+        degenerates = np.roots(trimmed)
+        all_roots[idx, : degenerates.size] = degenerates
+
+    return get_real_positive_smallest_roots(all_roots)

This preserves the fast path while restoring correct behaviour when the quartic term cancels.

📜 Review details

Configuration used: Path: .coderabbit.yaml

Review profile: CHILL

Plan: Pro

📥 Commits

Reviewing files that changed from the base of the PR and between f309cb8 and 83443bd.

📒 Files selected for processing (8)
  • pooltool/evolution/event_based/introspection.py (0 hunks)
  • pooltool/evolution/event_based/simulate.py (7 hunks)
  • pooltool/ptmath/roots/_quartic_numba.py (1 hunks)
  • pooltool/ptmath/roots/quartic.py (2 hunks)
  • tests/evolution/event_based/test_simulate.py (11 hunks)
  • tests/ptmath/roots/README.md (1 hunks)
  • tests/ptmath/roots/_generate_quartic_coeffs.py (1 hunks)
  • tests/ptmath/roots/test_quartic.py (1 hunks)
💤 Files with no reviewable changes (1)
  • pooltool/evolution/event_based/introspection.py
🧰 Additional context used
🪛 markdownlint-cli2 (0.18.1)
tests/ptmath/roots/README.md

19-19: Fenced code blocks should have a language specified

(MD040, fenced-code-language)

⏰ Context from checks skipped due to timeout of 90000ms. You can increase the timeout in your CodeRabbit configuration to a maximum of 15 minutes (900000ms). (4)
  • GitHub Check: test (3.13, 1.8.4, windows-latest)
  • GitHub Check: test (3.13, 1.8.4, ubuntu-latest)
  • GitHub Check: test (3.10, 1.8.4, windows-latest)
  • GitHub Check: test (3.10, 1.8.4, ubuntu-latest)

Comment thread pooltool/ptmath/roots/_quartic_numba.py Outdated
Comment on lines +425 to +431
if coeff[0] == 0.0:
return roots

a = coeff[1] / coeff[0]
b = coeff[2] / coeff[0]
c = coeff[3] / coeff[0]
d = coeff[4] / coeff[0]
Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🔴 Critical

Do not suppress lower-degree solutions when the quartic term vanishes

When coeff[0] == 0, solve immediately returns an array of zeros. Downstream, get_real_positive_smallest_roots treats that as “no positive root”, so any real solution of the remaining cubic/quadratic (e.g., grazing collisions where the t⁴ term cancels) gets lost and the simulation under-detects events. The legacy solvers handled this by falling back to lower-degree formulas; we need the same behaviour here instead of returning zeros.

Comment thread tests/ptmath/roots/README.md Outdated
Comment on lines +9 to +13
- `quartic_coeffs.npy` - Quartic polynomial coefficients (standard cases)
- `quartic_coeffcoeffs.roots.npy` - Ground truth roots for standard cases
- `hard_quartic_coeffs.npy` - Quartic polynomial coefficients (pathological cases)
- `hard_quartic_coeffs.roots.npy` - Ground truth roots for pathological cases
- `hard_quartic_coeffs.cases.npy` - Case labels indicating which pathological scenario each coefficient set represents
Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟡 Minor

Correct the hard coefficients file name

The second bullet currently reads quartic_coeffcoeffs.roots.npy, but the generated artifact (and tests) use quartic_coeffs.roots.npy. Please fix the typo so the README matches the actual files.

-- `quartic_coeffcoeffs.roots.npy` - Ground truth roots for standard cases
+- `quartic_coeffs.roots.npy` - Ground truth roots for standard cases
📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
- `quartic_coeffs.npy` - Quartic polynomial coefficients (standard cases)
- `quartic_coeffcoeffs.roots.npy` - Ground truth roots for standard cases
- `hard_quartic_coeffs.npy` - Quartic polynomial coefficients (pathological cases)
- `hard_quartic_coeffs.roots.npy` - Ground truth roots for pathological cases
- `hard_quartic_coeffs.cases.npy` - Case labels indicating which pathological scenario each coefficient set represents
- `quartic_coeffs.npy` - Quartic polynomial coefficients (standard cases)
- `quartic_coeffs.roots.npy` - Ground truth roots for standard cases
- `hard_quartic_coeffs.npy` - Quartic polynomial coefficients (pathological cases)
- `hard_quartic_coeffs.roots.npy` - Ground truth roots for pathological cases
- `hard_quartic_coeffs.cases.npy` - Case labels indicating which pathological scenario each coefficient set represents
🤖 Prompt for AI Agents
In tests/ptmath/roots/README.md around lines 9 to 13, there is a typo in the
second bullet: it reads `quartic_coeffcoeffs.roots.npy` but should match the
actual artifact `quartic_coeffs.roots.npy`; update that bullet to
`quartic_coeffs.roots.npy` so the README matches the generated files and tests.

@ekiefl ekiefl merged commit ab4b412 into main Nov 12, 2025
10 of 11 checks passed
@ekiefl ekiefl deleted the ek/1010-quartics branch November 12, 2025 03:52
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant