feat: add C++ backend for pdist_tangency_grad (batched pairwise tangency VJP)

[t-uda]

Closes #129

Summary

Adds a single-call batched C++ kernel pdist_tangency_grad to _tangency_cpp_impl.cpp that computes, over all C(N,2) pairs, the condensed tangency distances plus the per-pair gradient blocks dt_dp / dt_dq, using the same envelope-theorem + implicit-differentiation math as the Python reference (tangency_grad / solve_mu_gradients). It is exposed through the existing ctypes layer in _tangency_cpp.py and auto-dispatched from ellphi.grad.pdist_tangency_grad.

This follows the maintainer-approved answers to the issue's three questions:

1. Batched C++ entry point: yes — one call for all pairs, mirroring the forward pdist_tangency C export.
2. Naming / export surface: same as pdist_tangency — no new public API names and no separate build flag. pdist_tangency_grad(coefs) -> (dists, vjp) auto-selects the C++ kernel exactly when the C++ backend is available (same selection as the forward path), with an additional symbol probe (has_pdist_tangency_grad) so a stale cached library that predates the new export falls back to Python instead of failing.
3. API stability: the public signature and semantics are unchanged; the pure-Python implementation is kept as the reference fallback (ellphi.grad._pdist_tangency_grad_python), and the VJP accumulation (np.add.at scatter) stays in NumPy/Python. C++ returns distances + per-pair gradient blocks only.

Degenerate configurations (identical / concentric ellipsoids) keep raising ZeroDivisionError, matching the Python reference (new error-message mappings in _raise_backend_error).

Credit

The design and performance case come from the prototype by @koki3070 and collaborators (downstream project TDA-ML, local branch feature/pdist-tangency-grad-cpp, commit 9fcad0c). The prototype branch was not publicly accessible at implementation time, so this PR re-derives the kernel from the issue's description and the in-repo Python reference, matching this repository's C++/ctypes conventions. Thank you for the detailed issue and benchmarks!

Measured results (local, macOS arm64, N=120, d=2, 7140 pairs)

Path	forward + VJP
Python reference	~1191 ms
C++ kernel	~7.4 ms (~160x)

Numerical agreement vs the Python reference on random ellipse clouds (multiple sizes, 2-D and 3-D, including highly eccentric near-degenerate ellipses): distances agree to <= 1e-12 abs; gradients agree with max relative error ~1.8e-11 (max absolute error ~2.3e-9 at N=120, occurring only on entries of magnitude ~1e2; the test-suite tolerance assert_allclose(rtol=1e-10, atol=1e-10) passes on all covered clouds).

Timing is reported here only; tests do not hard-assert timing.

Implementation notes for reviewers

• The kernel reuses the existing helpers (solve_mu, pencil_into, cholesky_factor_into, solve_with_cholesky_into, quad_eval), so it works in both build modes: ELLPHI_USE_EIGEN=1 (verified: compiles and passes tests with cpp_linalg_kind='eigen') and the portable internal path (default, also verified).
• Per-pair math mirrors solve_mu_gradients step by step: pencil Cholesky, residual, dF/dmu = 2 r . K^{-1} r (with the same numpy.isclose threshold 1e-8), centre-Jacobian rows -K^{-1} rhs_k contracted against phi_x = -2 r (implemented as 2 (K^{-1} rhs_k) . r), then the envelope-theorem assembly with the 1/(2t) factor. Reviewers may want to double-check this contraction identity and the degenerate-error mappings.
• Rebuild trigger: the existing mechanism (source mtime + embedded TANGENCY_VERSION + linalg kind in scripts/build_tangency_cpp.py) already rebuilds the library on uv sync because the C++ source changed; no version-marker change is needed. For any environment that still loads an older same-version library, the new symbol probe falls back to Python gracefully (covered by tests).
• Minor behaviour-preserving refactor: the Python fallback now returns (dists, dt_dp, dt_dq) arrays and both paths share one vectorised np.add.at VJP (previously a per-pair Python loop closure).

Verification (all run in this branch)

• uv sync — OK (rebuilds backend)
• uv run python -m ellphi --build-info — OK (cpp_backend_available=True, checked with both cpp_linalg_kind='internal' and 'eigen')
• uv run black src tests scripts / uv run black --check src tests scripts — OK
• uv run flake8 src tests scripts — OK
• uv run mypy src tests — OK (no issues in 34 source files)
• MYPYPATH=src uv run stubtest ellphi --allowlist stubtest-allowlist.txt — OK
• uv run pytest — OK (259 passed; suite also re-run against the Eigen build)
• uv lock --check — OK
• uv run mkdocs build — OK (release-notes touch; pre-existing griffe warnings only)

🤖 Generated with Claude Code

https://claude.ai/code/session_01F1U8LhkDL8Ln2hgsMyFu2G

[codecov]

Codecov Report

❌ Patch coverage is 96.42857% with 2 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/ellphi/_tangency_cpp.py	93.10%	2 Missing ⚠️

📢 Thoughts on this report? Let us know!

[Copilot]

Pull request overview

This PR adds a batched C++ backend for ellphi.grad.pdist_tangency_grad, exposing it through the existing ctypes layer and auto-dispatching from the public Python API to dramatically speed up pairwise tangency VJP computation while preserving the existing public signature and Python fallback behavior.

Changes:

• Add a new C++ export pdist_tangency_grad that computes condensed distances plus per-pair gradient blocks (dt_dp, dt_dq) in one call.
• Update ellphi.grad.pdist_tangency_grad to dispatch to the C++ kernel when available (with a Python reference fallback) and to use vectorized np.add.at accumulation for the VJP.
• Add tests covering C++ vs Python numerical agreement, stale-library symbol probing, degenerate error cases, and fallback selection; update release notes and type stubs.

Reviewed changes

Copilot reviewed 6 out of 6 changed files in this pull request and generated 2 comments.

Show a summary per file

File	Description
tests/test_grad_cpp.py	Adds test coverage for the new C++ grad kernel, fallback behavior, and degenerate/stale-backend scenarios.
src/ellphi/grad.py	Adds Python reference helper and dispatches pdist_tangency_grad to C++ when available; refactors VJP accumulation to vectorized scatter.
src/ellphi/_tangency_cpp.pyi	Extends typing surface with has_pdist_tangency_grad and pdist_tangency_grad signatures.
src/ellphi/_tangency_cpp.py	Adds symbol probe and ctypes binding for pdist_tangency_grad, plus backend error mappings.
src/ellphi/_tangency_cpp_impl.cpp	Implements and exports the batched C++ kernel and per-pair gradient computation.
RELEASE_NOTES.md	Documents the new auto-dispatch behavior and performance improvement.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[t-uda]

LGTM

[t-uda]

Addressed both Copilot review findings in 6eab7c7:

1. VJP cotangent validation (src/ellphi/grad.py): the vjp closure now validates grad_dists up front. Chosen contract: the input must be broadcastable to shape (N*(N-1)//2,) — valid 1-D cotangents work as before and scalars broadcast to all pairs; anything else (wrong length, 2-D) raises a ValueError naming the expected shape. For the record, the pre-feat: add C++ backend for pdist_tangency_grad (batched pairwise tangency VJP) #134 loop implementation rejected scalars with low-level TypeError/IndexError, silently truncated over-long arrays, and accidentally accepted 2-D arrays with undocumented per-component semantics; the accidental acceptances are now explicit errors. Contract documented in the docstring, annotation widened to numpy.typing.ArrayLike, covered by tests (valid 1-D, scalar, wrong length, 2-D).
2. Stale-library guard (src/ellphi/_tangency_cpp.py): the internal binding now checks has_pdist_tangency_grad() and raises a RuntimeError with the module's standard rebuild guidance instead of a bare AttributeError; covered by a direct-call test.

Local checks re-run: black / flake8 / mypy / stubtest / pytest (264 passed) all green.

🤖 Generated with Claude Code

[koki3070]

Downstream benchmark (TDA-ML prototype vs this PR)

We compared our local prototype (feature/pdist-tangency-grad-cpp, 9fcad0c) against this PR branch (6eab7c7) on the same fixed inputs in TDA-ML.

Environment: Linux x86_64 (AMD Ryzen Threadripper PRO 7995WX), Python 3.12.13, warmup=3 / repeats=10 (median reported).

Speed (C++ kernel, N=120, d=2, 7140 pairs)

	Prototype (9fcad0c)	PR (6eab7c7)	Python reference
C++ kernel	8.84 ms	9.17 ms	—
forward + VJP (public API)	9.49 ms	9.81 ms	~1.8 s
speedup vs Python	~200×	~196×	—

Also checked TDA-ML elongate workload (MNIST + local-PCA teacher, N=100): C++ 6.4–7.0 ms vs Python ~1.2 s (~175–193×).

Accuracy (same coefs for both builds)

• 2-D cases + TDA-ML elongate coefs: C++ outputs match bitwise between prototype and PR (max abs error 0 on dists, dt_dp, dt_dq).
• 3-D N=50: max abs error ≤ 1.3×10⁻¹¹ on dt_dq (prototype built with eigen, PR with internal; negligible).
• Both vs Python reference (N=120): distances exact; gradient max relative error ≤ 1.2×10⁻¹².

Conclusion

From a downstream perspective, this PR’s re-implementation matches our prototype numerically and delivers the same order of speedup on Linux (~200× at N=120). Looks good to merge from our side — thank you for the clean upstream port!

Additional cases (timing median ms)

case	Python (proto)	C++ (proto)	C++ (PR)
random_n10_d2	10.9	0.06	0.07
random_n30_d2	106.2	0.55	0.61
random_n50_d3	344.5	2.02	2.27
eccentric_n30_d2	106.9	0.90	0.79
tda_ml_elongate_n100	1222.9	6.96	6.38

[t-uda]

@koki3070 ありがとう．random ケースだけ proto 版の方がほんの若干だけ速いみたいですが，ほかはそうでもないところだけ見るとこれはただの誤差ですかね？ 今の PR 実装にだけ何かしらかのオーバーヘッドがあったりしないかだけ少し気になりました．

オーバーヘッドがなにかあったとしても影響は軽微そうなので，様子を見てからマージします．