Herget IOD

[awilson110]

Fixes #311

Implemented the Herget IOD method, using assist variational particles for calculating partial differentials.

Review Checklist for Source Code Changes

• Does pip install still work?
• Have you written a unit test for any new functions?
• Do all the units tests run successfully?
• Does Layup run successfully on a test set of input files/databases?
• Have you used black on the files you have updated to confirm python programming style guide enforcement?

[matthewholman]

Four things before merge. Line numbers are in src/layup/utilities/herget_iod.py unless stated.

1. :65 mutates the caller's Observation objects and never restores the epochs -- obs holds references, so the LM fit afterwards runs on times shifted by ~0.23 d at the 40 au guess.
2. :253 the inner while has no iteration cap, so non-convergence hangs rather than fails (cf. Unusually long orbitfit execution time for certain MBAs #465).
3. :88 flag = 0 is set even when the outer loop exits on max_iterations.
4. :91 cov = [0.01]*36 is singular -- every entry equal gives one nonzero eigenvalue.

Smaller: live print at :77 and test_herget_iod.py:127; GMtotal at :291 duplicates layup.constants.MU_SUN; :172 discards the find_velocity result for the rn endpoint; docstrings still list args/aux.

[matthewholman]

Correction to item 1 above. I asserted the fit afterwards is corrupted; I had
read that, not run it. Having now measured it, the second half is wrong.

The mutation is real: 61 of 61 epochs shifted by -5.44 h (-0.2266 d, the light
time at the converged ~39 au) and not restored.

But it does not degrade the fit. Against the JPL state in your own
test_whole_herget_method:

iod	flag	|dr| from JPL (au)
gauss	0	0.0864
herget	0	0.0014

Herget lands 64x closer than the Gauss-seeded fit, so nothing is being spoiled
here.

So item 1 stands only as a side effect on caller-owned state: anything reusing
those observations afterwards silently gets emission times rather than reception
times. Restoring the epochs before returning is still worth doing, but it is
hygiene, not a blocker. Items 2-4 are unaffected.

Incidentally the gauss error above exceeds the atol=1e-2 that test asserts.

[matthewholman]

Ran this on 84 real Rubin-only arcs (MPC obscode X05, median 15 d, ~12 obs).

The initial range is the main issue. rho_1 = rho_n = 40 is a TNO value.
Every arc inside 5 au failed and every one outside converged. Setting it to
2.0 au takes convergence from 9.5% to 77%, and the distant objects still
converge -- so it is not a trade.

Where both converge, herget and gauss agree on a to a median |da/a| of 1.7e-6.
The method is fine.

The remaining ~23% is a degeneracy, not a bug. a1 and a2 come out nearly
antiparallel (corr -0.998 typically, reaching 1-|corr| ~ 1e-16 on the failures),
the normal matrix collapses, and the unregularised step runs to ~1e11 au until
the epochs leave the ephemeris. That is the short-arc range degeneracy
(Milani et al. 2004, CMDA 90, 57) -- it wants damping or a step limit.

One more: :134 uses a 1 au variational displacement, which gives a partial
~40% low at 1 au. Scale it to the range -- and the find_velocity tolerance
with it, or the difference drops below the solver's own noise.

[matthewholman]

Where this matters most: the arcs the current default cannot fit at all.

I took the 55 Rubin-only arcs (X05) whose iod=gauss fit exceeded 420 s and
returned nothing, and re-ran them:

iod	converged	median time
herget (rho0 = 2 au)	41/55 (75%)	0.15 s
auto	0/55	all timed out at 420 s

Fitted semimajor axis: median 2.75 au, range 1.37-3.23 -- main belt.

auto does not help because the BK-IOD fallback only fires when the Gauss roots
fail; here a root is found and then grinds on a close-Earth passage (#465), so
the fallback never runs.

So this is not a redundant third IOD -- it covers the regime where the default
fails hardest, and it does so ~3000x faster than the failure costs. Relevant to
#486 as well: the ladder probably wants herget in it, tried early when the sky
motion is fast, since that predicts this failure before any fitting.

Caveat: this used rho_1 = rho_n = 2.0, not the 40 in the branch. With 40 none
of these converge. Worth reproducing independently before relying on it.

[matthewholman]

On the degeneracy above -- a concrete suggestion.

Two lines would convert the runaway into an honest non-convergence:

1. Levenberg damping on the 2x2 normal matrix, S11*(1+lam) and
   S22*(1+lam) in the determinant, with lam raised when a step increases
   the residual and lowered when it decreases. Same scheme the main fitter
   already uses.
2. A step cap, |delta_rho| <= 0.5 * rho. Crude, but on its own it would
   have caught every failure I saw -- nothing would reach 1e11 au.

Worth knowing the shape of it: since a1 and a2 are nearly antiparallel, the
badly determined direction is rho_1 + rho_n, i.e. the mean range, while
rho_1 - rho_n (the range rate) is fine. Damping refuses to move along the bad
direction softly; dropping the small eigenvalue refuses completely; ranging
methods refuse and enumerate it instead. Same recognition, increasing force.

Damping plus the cap is the minimal fix and does not change what the method
claims to determine.

[mschwamb]

@matthewholman maybe this is your text in the comments here and not copy and pasting output from Claude but I can't fully parse what the content of all the comments here.

What are the 55/84 Rubin orbits referred to in these comments @awilson110 can replicate and learn?

[awilson110]

@matthewholman I've made a fix for points 1-3 (and the smaller points too), but I'm stuck on both point 4 and the degeneracy suggestion.

For point 4: the code only calculates one eigenstate (hence cov being singular is no problem), should I change this? How would I calculate for multiple eigenstates? The only way I could see is to run the code for a range of different initial guesses and see if they converge differently, is that what I should do?

For the degeneracy: I'm not familiar with the notation, what are S11*(1+lam), S22*(1+lam), a1, a2?

[matthewholman]

Adam — drop point 4 and the damping suggestion. Neither is where Herget earns its place. (I can explain the notation another time.)

It turns out that gauss actually covers a much larger fraction of the cases, but Herget is much faster in some cases when gauss leads to a very slow result.

So the thing to work on is the initial range. rho_1 = rho_n = 40 is hardcoded, and the 2 au that works on other objects was tuned on the main belt. Can you make it geometry-derived, or a short ladder, and test on a near-Earth and a main-belt object as well as the TNO?