Skip to content

Fix quadrature noise in linking number diagnostic - #63

Open
rogeriojorge wants to merge 3 commits into
mainfrom
fix/integer-linking-number
Open

Fix quadrature noise in linking number diagnostic#63
rogeriojorge wants to merge 3 commits into
mainfrom
fix/integer-linking-number

Conversation

@rogeriojorge

@rogeriojorge rogeriojorge commented Aug 24, 2026

Copy link
Copy Markdown
Member

Summary

  • classify each coil pair by rounding its absolute Gauss integral to the integer linking number, matching SIMSOPT's compute_linking_number
  • explicitly stop gradients through this topological diagnostic so quadrature noise cannot exert a spurious optimization force
  • expose focused internal helpers for the signed per-pair quadrature values and per-pair integer classifications; the public API still returns their sum
  • use the exact |r1-r2|^3 denominator used by SIMSOPT, while safely masking padded block entries
  • test known topologies, every symmetry-expanded coil pair, blockwise evaluation, curve conventions, and derivative behavior against both an independent NumPy reference and SIMSOPT itself

Motivation

For disjoint closed curves, the Gauss linking number is integer-valued. The previous implementation returned the raw quadrature approximation, so unlinked circular coils could report a small nonzero loss. More seriously, differentiation of that residual could produce a nonzero gradient and perturb the coils even though linking number cannot vary continuously without an intersection.

SIMSOPT handles this by rounding each pair's absolute Gauss integral before summing and provides no derivative for this topological diagnostic. This PR adopts the same semantics in JAX using per-pair round and stop_gradient.

The former 1e-12 denominator offset was also removed because it changes the Gauss integral and prevented a point-for-point comparison with SIMSOPT. Invalid padded entries are now assigned a safe denominator and explicitly zeroed. A true zero separation remains singular, which is appropriate because linking number is undefined when curves intersect.

Test coverage

The non-optional tests contain a NumPy translation of SIMSOPT's C++ loop, so CI verifies the same algorithm without adding SIMSOPT as a dependency. They check:

  • all 66 pairs of 12 nfp=3, stellsym=True symmetry-expanded initial circular coils
  • unlinked circles classify as zero per pair and in aggregate despite raw quadrature residue
  • a Hopf link classifies as one
  • signed-integral behavior under orientation reversal and coil-order exchange
  • block_size=None and non-divisible block_size=16 with 65/97 quadrature points
  • exactly zero JAX gradient

The optional integration tests use pytest.importorskip("simsopt"). When SIMSOPT is available they additionally verify:

  • ESSOS and SIMSOPT gamma and gammadash arrays point-by-point for identical Fourier dofs, including every symmetry copy
  • raw per-pair ESSOS integrals against the NumPy/C++ reference
  • integer classifications against a SIMSOPT LinkingNumber for every pair
  • aggregate totals against SIMSOPT
  • ESSOS's exact-zero gradient alongside SIMSOPT's zero derivative

Local results

Using the SIMSOPT checkout installed from /Users/rogeriojorge/local/simsopt_test:

  • python -m pytest tests/test_linking_number_reference.py tests/test_linking_number_simsopt.py -q — 8 passed
  • python -m pytest tests/test_objective_functions.py tests/test_linking_number_reference.py -q — 17 passed
  • python -m pytest -q — 90 passed, 1 xfailed

No production dependency was added. On CI without SIMSOPT, only the direct integration module is skipped; the independent NumPy reference coverage remains mandatory.

@codecov

codecov Bot commented Aug 24, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.

Files with missing lines Coverage Δ
essos/objective_functions.py 99.08% <100.00%> (+0.34%) ⬆️

... and 1 file with indirect coverage changes

@rogeriojorge

Copy link
Copy Markdown
Member Author

Added SIMSOPT-compatible downsample support in 0fd8379. The integration step is scaled by the factor, invalid/non-divisible values are rejected, and the existing blockwise path operates on the reduced point set. Tests compare raw per-pair integrals and exact integer classifications against both the NumPy translation of SIMSOPT C++ and installed SIMSOPT for factors 1, 2, 3, 4, 5, and 6, covering expanded unlinked circles and a Hopf link. Full local suite: 99 passed, 1 xfailed. Benchmark for 56 expanded coils / 80 source points: 6.61 ms at downsample=1, 1.33 ms at 2, and 0.45 ms at 8; all report linking number 0.

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