Skip to content

P2.1 Spherical potentials: backend-agnostic namespace-swap#893

Merged
jobovy merged 4 commits into
feat/backendsfrom
feat/P2.1-spherical
Jun 6, 2026
Merged

P2.1 Spherical potentials: backend-agnostic namespace-swap#893
jobovy merged 4 commits into
feat/backendsfrom
feat/P2.1-spherical

Conversation

@jobovy

@jobovy jobovy commented Jun 4, 2026

Copy link
Copy Markdown
Owner

Part of the P2.x potential namespace-swap series (→ #892 integration branch).

Migrates the spherical family to backend-agnostic compute (numpy/jax/torch) via the PR0 convention (xp = get_namespace(...); private compute layer only).

Migrated: SphericalPotential (the 3D→1D reduction layer — covers the whole family), Dehnen(Core)Spherical, Hernquist, Jaffe, NFW, TwoPowerSpherical (base dens/z2deriv), PowerSpherical/Kepler, HomogeneousSphere, SphericalShell, Burkert.
Deferred to Pspecial (scipy.special/interpolate): TwoPower _evaluate (hyp2f1/gamma), NFW _evaluate (xlogy), the _surfdens of Hernquist/Jaffe/NFW/Burkert (data-dependent complex-sqrt branch), Einasto, interpSpherical.
Hazards fixed: torch lacks power**; in-place out[r2==0]=-inf→safe-denominator xp.where; piecewise if r<a→shape-polymorphic xp.where; numpy.fabsxp.abs (a real cross-namespace bug caught by the test).
Tests: tests/test_backend_spherical.py — 95 passed, 15 skipped (numpy/jax/torch parity rtol 1e-12 + grad-vs-FD + grad==−Rforce). Numpy path byte-identical (max diff 0.0 on spot checks).

🤖 Generated with Claude Code

@codecov

codecov Bot commented Jun 4, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 99.92%. Comparing base (141ebfe) to head (9ac1806).
⚠️ Report is 1 commits behind head on feat/backends.

Additional details and impacted files
@@               Coverage Diff               @@
##           feat/backends     #893    +/-   ##
===============================================
  Coverage          99.92%   99.92%            
===============================================
  Files                228      228            
  Lines              34213    34461   +248     
  Branches             712      716     +4     
===============================================
+ Hits               34186    34434   +248     
  Misses                27       27            

☔ View full report in Codecov by Harness.
📢 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.

@jobovy jobovy force-pushed the feat/P2.1-spherical branch from 7628b5a to 9432cca Compare June 5, 2026 16:54
jobovy and others added 2 commits June 5, 2026 20:46
The P2.1 migration of PowerSphericalPotential._evaluate switched to an
xp.where-based form, so _evaluatePotentials now returns a 0-d numpy
scalar (with .shape == ()) instead of a Python float for scalar input.
This made fE/fQ take the array (reshape) branch instead of the scalar
`return out[0]` branch, dropping coverage of that line in
isotropicPowerLawdf.py (L126) and osipkovmerrittPowerLawdf.py (L152)
and failing codecov/project (+2 misses).

Add direct scalar-input assertions (matching the existing
constantbetaPowerLawdf pattern) to restore 100% coverage. Test-only
change; numpy code path is byte-identical.

Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com>
The backend rewrite added safe_h = where(h==0, 1, h) to avoid 1/0, but at the
exact shell edge R==a the projected surface density genuinely DIVERGES (the
1/sqrt(a^2-R^2) limb singularity), so 1/(2 pi a h) -> inf is the physically
correct value and what the original numpy code returned. Masking h there
silently replaced the divergence with a finite ~1/(2 pi a). Drop safe_h and use
h directly: inf at R==a (matches original, byte-identical for all R), the only
consequence being a singular gradient exactly at the edge (which is real). The
R>a branch remains NaN-safe via safe_arg (sqrt(negative) guard) above.

Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com>
@jobovy jobovy marked this pull request as ready for review June 6, 2026 13:29
Migrate the remaining numpy-only / value-branching compute methods of the
spherical-potential family to the galpy.backend (xp) namespace layer, so they
run and differentiate under numpy / jax / torch:

- BurkertPotential._surfdens, ._revaluate
- HernquistPotential._surfdens, JaffePotential._surfdens, NFWPotential._surfdens
- NFWPotential._evaluate, NFWPotential._mass
- SphericalPotential._mass
- PowerSphericalPotentialwCutoff._dens, ._ddensdr, ._d2densdr2

The four cuspy _surfdens methods used numpy.sqrt(R^2-a^2+0j) complex arithmetic
plus an `if Rma == 0` scalar branch for the R == a removable singularity. These
are rewritten with a complex sqrt built from a NaN-safe real argument
(xp.astype(safe, complex128)), xp.real for the generic branch, and an xp.where
selecting the closed-form edge limit at R == a. Per the safe-denominator rule
(jax/torch where evaluates both branches), every 1/Rma, 1/(R^2-a^2),
1/(r^2-a^2) in the dead branch is guarded so it cannot inject NaNs into
reverse-mode autodiff.

NFW._evaluate and Burkert._revaluate used scipy.special.xlogy purely for its
0*inf -> 0 convention as r -> infty (the prefactor 1/r, 2/x is never zero for
finite r); these are reproduced exactly with a NaN-safe xp.where, so no scipy
dependency remains. All migrated numpy paths are bit-identical (atol=0) to the
originals, including the R == a edge and the r == 0 / r == inf limits.

Left numpy-only (genuine scipy.special, no backend-agnostic replacement yet):
base TwoPowerSphericalPotential._evaluate/_Rforce/_zforce/_R2deriv/_Rzderiv/
_mass (hyp2f1, gamma), PowerSphericalPotential._surfdens (hyp2f1),
PowerSphericalPotentialwCutoff._evaluate/_Rforce/_zforce/_R2deriv/_z2deriv/
_Rzderiv/_mass (gamma/gammainc/hyp1f1).

tests: extend test_backend_spherical.py with numpy-vs-jax-vs-torch parity for
the migrated _surfdens (grid spanning R<a, R==a, R>a), an edge-finiteness
check, grad-vs-finite-difference off the edge plus grad-finiteness at the edge
(exercising the safe-where guards), the PowerSphericalPotentialwCutoff closed-
form densities, and add _evaluate to the NFW and Burkert AD cases.

Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com>
@jobovy jobovy merged commit 562fe31 into feat/backends Jun 6, 2026
108 checks passed
@jobovy jobovy deleted the feat/P2.1-spherical branch June 6, 2026 20:13
jobovy added a commit that referenced this pull request Jun 10, 2026
Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com>
jobovy added a commit that referenced this pull request Jun 12, 2026
Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com>
jobovy added a commit that referenced this pull request Jun 13, 2026
Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com>
jobovy added a commit that referenced this pull request Jun 15, 2026
Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com>
jobovy added a commit that referenced this pull request Jun 17, 2026
Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com>
jobovy added a commit that referenced this pull request Jun 19, 2026
Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com>
jobovy added a commit that referenced this pull request Jun 19, 2026
Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com>
jobovy added a commit that referenced this pull request Jun 25, 2026
Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com>
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