orbit: the V_c/V_e eccentricity parameterization (fitvcve)

[jdeast]

fitvcve: true now works, per orbit: it samples V_c/V_e -- the ratio of the
speed a planet would have on a circular orbit of the same period to its speed at
transit -- and the direction of omega, instead of sqrt(e)cos(omega) /
sqrt(e)sin(omega). The reason is Eastman (2024, PASP; arXiv:2309.14410): over 330
simulated systems, a transit-only fit in the sqrt(e) pair recovers eccentricities
that are measurably wrong, while the transit duration constrains
V_c/V_e = sqrt(1 - e^2)/(1 + e sin omega) tightly. Review 8.8.3; the stubs it
replaces raised.

What the paper does not address, and this does

The paper picks a root of the quadratic inversion with a discrete sign parameter
S, and rejects combinations with no real root. Both are fine for DE-MC and
useless to a gradient sampler -- S is piecewise constant with a logp that jumps
at S = 1, and a rejection is a wall with nothing to follow back out. So:

• the likelihood is marginalized over both roots (System._add_branch_mixtures:
  one logsumexp over the 2^k root combinations of k V_c/V_e orbits, logged,
  warning past two). It substitutes into the RV-level term sum and adds
  logsumexp(...) - L_ref, which cancels the copy PyMC already added -- so no
  data component knows this is happening, a term independent of ecc factors out
  identically, and one that depends on it (the Jacobian, the collision bound, a
  future orbit-crossing penalty) is replicated per branch, making each branch's
  weight log w + log|J| + barriers: the form review 8.4.4 specifies for folded
  likelihoods, for free;
• the discriminant is floored inside the sqrt (hard: a NaN is unbuildable)
  while the unfloored value drives a soft_lower_bound (soft: the floored one
  is flat across the whole imaginary region, and a flat penalty has no gradient).

Also here: role 3's deferred build pass (the last piece of 8.2.1, which 8.8.3
owned), vcve at rank: 5 so the relaxation engine rewrites it rather than
pushing back on ecc/omega, and parameterization.merge_options.

Five bugs the verification caught, all mine

1. A dropped factor of 2 in the discriminant. B^2 - 4AC = 4(1 - x^2 cos^2 w),
   so the roots need (-B +/- 2 sqrt(...))/2A. The first draft's roots were real,
   in range, and wrong -- caught by round-tripping the inversion against eq 4 over
   an (e, omega) grid, which is now a test for exactly that reason.
2. "Primary = the lower root" is wrong. The roots' product is C/A, so below
   V_c/V_e = 1 they have opposite signs and the physical one is always the
   upper. Preferring the lower would have pinned the forward model at e = 0,
   with no gradient, across that whole region.
3. orbit.tp consumed a reported element. It reads secosw/sesinw, which a
   V_c/V_e orbit reports rather than samples, so it would have read the
   pre-patch placeholder: a silently wrong periastron. Fixed with per-orbit
   (e, omega) expressions for tp, esinw and ecosw -- and
   calc_tp_from_ecc is strictly better behaved besides (writing the true
   anomaly at conjunction as f = pi/2 - omega cancels the sqrt(e) factor,
   removing both of calc_tp's e = 0 singularities, review 1.8.2). This is why
   System._validate_reported_not_consumed now exists; it found esinw/ecosw
   the moment it caught tp.
4. The Jacobian had the wrong sign. V_c/V_e is the sampled coordinate, so
   p(e) ~ |d(V_c/V_e)/de| -- the paper's high-eccentricity bias -- and
   flattening it means adding the reciprocal. Adding the derivative doubles
   the bias, and the finite-difference test passes either way, so the direction is
   now pinned by measuring the implied density on e for flatness plus a
   call-site test against that. (EXOFASTv2's convention is deliberately not
   restated in the comments: it could not be checked from this tree.)
5. vcve's upper bound of 2.0 capped e at 0.6 for an orbit transiting near
   apastron (eq 4 at omega = -90 deg is sqrt((1+e)/(1-e))) -- silently
   excluding exactly the eccentric transiting planets the parameterization exists
   to measure. It is now the image of the eccentricity cap,
   sqrt((1+MAX_ECC)/(1-MAX_ECC)) = 141.4.

A sixth was in the mixture itself: with two V_c/V_e orbits, both branches replace
the same ecc node, and merging their replacements into one dict kept only the
last -- silently marginalizing over 3 of the 4 combinations. They are applied in
sequence now (which composes, because each replacement reads the node it
replaces), checked at declaration, and measured against a closed-form four-branch
marginalization whose branch offsets are deliberately small: with a far-away
branch the logsumexp is dominated by one term and the missing combination is
invisible.

Still guarded

The chord half (fitchord: true, cosi's from_b, the orbit-level b),
whose physics is undefined. A fitchord that merely follows fitvcve does not
raise -- it was not asked for. The transit-only default (both halves on for an
orbit constrained by transits alone) flips when that lands, so the pair the paper
validated turns on together and the shipped transit-only examples move exactly
once.

Verification

• All 19 shipped example configs keep their exact start logp, free-RV set and
  sampled structure. The only difference is the new orbit.vcve row, now
  REPORTED on every orbit, and its values sanity-check the forward relation
  (hd80606's e ~ 0.93 gives 1.82, hat3's near-circular orbit 0.99986,
  GaiaBH1's e ~ 0.45 gives 0.81).
• Full suite green (2435 passed). tests/test_vcve.py is 79 tests; the fitvcve
  guards in tests/test_wip_guards.py are rewritten for the new behavior.
• A fitvcve transit fit really samples under nuts_sampler="numpyro" -- three
  separate places in this feature could produce a NaN only in a JAX backward
  pass.

🤖 Generated with Claude Code