Microlensing start values invert an underdetermined system: mu_rel/pi_E direction is picked arbitrarily

[jdeast]

Split out from #92, which made this arbitrary choice deterministic without making it right.

The relations are correct; the inverse problem is not determined

The forward physics in components/mulensing/symbolic_physics.py is right:

relation	check
mu_rel_mag**2 = mu_ra_rel**2 + mu_dec_rel**2	Pythagorean
pi_E_N = (pi_rel/theta_E)*(mu_dec_rel/mu_rel_mag)	pi_E || mu_rel; North<->dec
pi_E_E = (pi_rel/theta_E)*(mu_ra_rel/mu_rel_mag)	East<->RA
pi_rel = KAPPA*m*(pi_E_N**2 + pi_E_E**2)	from theta_E = sqrt(kappaMpi_rel), pi_E = pi_rel/theta_E

Forward (physical -> observable) is single-valued and fine. Initialization runs it backwards, and there the system is underdetermined.

What the seed actually contains

For examples/DC2018_128, mmexofast.json carries only:

parameters: t_0, u_0, t_E, s, alpha, rho, q
sigmas:     t_0, u_0, t_E, log_rho, log_s, log_q, alpha

No pi_E, in either fit entry. So with t_E seeded and theta_E from mass/distance, the engine knows only the magnitude mu_rel_mag = theta_E / t_E. Splitting that into (mu_ra_rel, mu_dec_rel) is one equation in two unknowns -- a circle of solutions. The pi_E vector direction is undetermined for the same reason.

What the engine does with that

It assigns the entire magnitude to whichever component it reaches first and leaves the other at (numerically) zero. Observed start values at the pinned GOOD_RAW:

• radish: mu_ra_rel = -11.4113305, mu_dec_rel = 0.0, pi_E_N = 0.0, pi_E_E = -0.22019125
• Ubuntu 26.04 (pre-config: sort component discovery so start values cannot depend on the filesystem #92): the transpose, with mu_ra_rel = 4.7e-125

Both satisfy the magnitude equation exactly. They are different physical models.

Why this matters even now that it is deterministic

1. The direction is observable. With annual parallax modeled, the sky direction of mu_rel affects the light curve -- which is precisely why the two branches gave chi2/N 9.96 vs 268.93 on the same data.
2. Even the good branch is a poor start. chi2/N ~ 10 at a seed meant to sit near the MMEXOFAST solution suggests the arbitrary direction is leaving real information on the table on every machine, not just the one that exposed the bug.
3. Exactly zero is a bad place to start. Putting the whole magnitude in one component makes the other exactly 0, so pi_E_N = 0 begins at a cusp of the normalization rather than in the interior.

Options

1. Seed the direction from the galactic model, which already computes expected proper motions (components/galacticmodel/galacticmodel.py:364-365) -- a physically motivated split instead of "all in RA".
2. Stop solving for it. Treat the mu_rel position angle as a free parameter with a prior, since the seed genuinely does not determine it.
3. Have MMEXOFAST fit parallax and seed pi_E_N/pi_E_E for events with parallax signal.

(1) looks cheapest and strictly better than the status quo; (2) is the honest representation of what is known.

Unrelated thing noticed nearby

galacticmodel.py:364 passes pm_ra_cosdec=pm_ra, i.e. it reads star.pm_ra as mu_alpha*cos(delta), while mu_ra_rel = lens_pm_ra - source_pm_ra is convention-agnostic. Self-consistent only if every consumer of star.pm_ra uses the cos(delta)-multiplied convention. Not audited -- flagging, not claiming a bug.

Reproducing

diag_mulens.py (in #92's discussion) prints versions, all resolved initvals/scales, every parameter's physical value at GOOD_RAW, and the mulens residual summary. Run it on two machines and diff.

🤖 Generated with Claude Code

[jdeast]

Measured: the family is NOT degenerate, and a correction to the framing above

Raised in review: if the system is underdetermined, shouldn't the whole family give the same light curve? That was the rationale for picking arbitrary values. Tested directly.

Built DC2018_128 twice with the lens proper motion rotated 90 degrees -- exactly the transpose the ordering bug produced -- and evaluated both at the seed (raw = 0), so no pinned raw offsets confound it:

	A: hint on pm_ra	B: rotated 90 deg
|mu_rel| heliocentric	11.411331	11.411331
|mu_rel| geocentric	11.489481	11.309388
t_E	18.047 d	18.334 d
pi_E (N, E)	(-0.220183, +0.001959)	(-0.001513, -0.220186)
chi2 (870 pts)	1316.05	1239.21

The magnitudes are identical to seven digits, so it genuinely is one family. But the light curve distinguishes its members, for two reasons:

1. The helio -> geo conversion is not rotationally symmetric. mu_rel_geo = mu_rel_helio - (Earth transverse velocity term) is a vector subtraction, so rotating mu_rel_helio changes |mu_rel_geo| (11.4895 vs 11.3094) and therefore t_E (18.05 vs 18.33 d). t_E sets the event timescale directly.
2. pi_E is parallel to mu_rel, so the rotation rotates the parallax vector. With annual parallax modeled, that direction is observable.

Net Delta chi2 ~ 77 over 870 points at the seed. Picking an arbitrary member is therefore not safe.

Correction to the numbers in the issue body

The -945.57 -> -113614.65 blow-up is not this effect. At the seed the two arrangements differ by only Delta chi2 ~ 77. The catastrophic gap appears at the pinned GOOD_RAW because the same raw offsets are applied to transposed initvals, so the physical vectors stop being a pure rotation -- |mu_rel| became 18.27 vs 28.88. The regression test amplified a real but modest seeding difference.

So "two different physical models from one config" overstates it. The accurate statement: nondeterministic seeding produced starts differing by Delta chi2 ~ 77, which the pinned-raw-point test amplified into a 225,000 chi2 gap. #92 is still a real bug -- a start that depends on the filesystem is indefensible regardless of magnitude -- but the severity framing there should be read with this in mind.

Direction chosen

Option 1 (seed from the galactic model's expected proper motion) is preferred over "pick the lowest-chi2 member". Lowest-chi2 works but is a mini-optimization at seed time that only searches the circle as parameterized. The galactic prior is physically motivated, is a genuine prior draw rather than an arbitrary corner, and is a better MCMC start.

It also fixes something both current branches share: each leaves one component at ~0 (pi_E_E = +0.002, pi_E_N = -0.0015), so a parameter starts pinned at zero, where pi_E_N sits at a cusp of the normalization.

Reproducer: scripts/diag_mulens.py for the cross-machine diff; the rotation test above is a two-build variant holding |mu_rel| fixed and evaluating at raw = 0.

[jdeast]

MMEXOFAST fits with pi_E = 0, which constrains what "correct seeding" can mean

Important context raised in review, and it changes how the numbers here should be read.

MMEXOFAST assumes pi_E = 0 when it produces the best fit we seed from. So the seeded t_0, u_0, t_E, s, alpha, rho, q describe a no-parallax model. This model then derives pi_E = pi_rel / theta_E from the physical parameters and gets something of order 0.2 -- not small. The seed is therefore internally inconsistent by construction: light-curve parameters fitted without parallax, evaluated in a model that applies a substantial one.

Consequences:

1. chi2/N > 1 at the seed is expected, not a bug. The measured 1.42-1.51 at raw = 0 is partly this mismatch. An earlier comment in this thread cited "even the good branch starts at chi2/N ~ 10" as evidence the seeding wastes information; that figure was at the pinned GOOD_RAW, and the seed itself is 1.4-1.5. Neither number should be read as a target of ~1, because we are not fitting MMEXOFAST's model.
2. We cannot reproduce MMEXOFAST's fit exactly, and should not try. Its chi2 is achievable only with pi_E = 0; matching it would mean giving up parallax, which is one of the things this model exists to measure.
3. It sharpens the choice of seeding target. Options for the direction, given the magnitude |pi_E| = pi_rel/theta_E is already fixed by the distances and mass:
   • Seed the direction from the galactic model's expected proper motion (the agreed direction). Physically motivated, avoids the exactly-zero component, good MCMC start.
   • Or start pi_E near zero, consistent with the fit the seed came from, and let the sampler find the parallax. Cheap, and self-consistent with t_E as seeded -- but it starts at the cusp / edge case this issue already flags as bad, and abandons the prior information.

The first still looks right. The second is worth noting only because it is the option that keeps t_E self-consistent with its own source, and the tension between "self-consistent with the seed" and "consistent with the prior" is unavoidable while the seed comes from a pi_E = 0 fit.

On the degeneracy measurement

scripts/diag_rotation.py (added in #92) now pins all four proper-motion components rather than just the lens's. Pinning only one lets the engine re-derive the other, and |mu_rel| drifts 11.411 -> 11.616 -- a 1.8% contamination that would quietly invalidate the comparison. With correct pinning the numbers reported earlier hold: |mu_rel| identical to seven digits, chi2 1239.21 vs 1316.05, delta 76.84. The script asserts the magnitudes match and aborts rather than reporting a contaminated delta.

[jdeast]

Reopened

This was auto-closed when #99 merged: that PR's body opened with "Closes #93 for the case it was filed about", and GitHub acts on the keyword and ignores the qualifier. The general problem is not solved. Shelving deliberately for now, with the state written down.

What #99 actually settled

star.pm_ra/star.pm_dec are now seeded at the galactic model's kinematic-prior mean (galacticmodel/physics.py:expected_proper_motion) instead of being left at defaults for the t_E constraint to absorb one component of. On DC2018_128 that takes |mu_rel| from 11.41 to 4.88 mas/yr -- the canonical microlensing value -- with t_E unchanged at its seed and theta_E/lens mass absorbing the difference, which is the standard microlensing chain and needed no special-casing.

The MMEXOFAST automated path is covered. That is the case this issue was filed about: MMEXOFAST fits with pi_E = 0 and emits no pi_E, so nothing in its output says which way mu_rel points. Verified live -- examples/DC2018_128 logs Multi-seed sampling: solved 2 seed start points and then [lens] star.0 proper motion seeded at the thin_disk prior mean for 4000 pc: pm_ra=-0.346, pm_dec=-1.323 mas/yr, and examples/DC2018/run_event.py:180 deliberately omits the microlensing start values so auto-init triggers.

What is still open

Seeding is gated off where a config already implies the proper motion, so on those events the engine still picks the direction itself -- deterministic since #92, but still arbitrary:

example	why gated	status
ob140939	pi_E_N/pi_E_E measured (Yee+2015)	direction data-determined; prior skipped
ob161003	two sources	direction still invented by the engine

The gate exists because seeding unconditionally made both of these much worse at the seed -- chi2/N 3.04 -> 179.11 and 1.72 -> 3.91. Root cause: hinting both pm components imposes a magnitude as well as a direction, and where the data already fix |mu_rel| (via t_E plus theta_E, or a measured pi_E) the prior fights them and something measured bends.

Why the obvious fix is not available

Filling only the direction and leaving the magnitude to the data would serve every case at once. Provenance cannot express that: rank is per-symbol and one-dimensional, and the direction is not a symbol. Gating on probe_derivable only patches ob140939, not ob161003.

Measured for the record, holding |mu_rel| fixed at 11.411331 and rotating 90 degrees (scripts/diag_rotation.py, at raw = 0):

	|mu_rel| helio	|mu_rel| geo	t_E	chi2
mu_rel in RA	11.411331	11.309388	18.334 d	1239.21
rotated 90 deg	11.411331	11.489481	18.047 d	1316.05

delta chi2 = 76.84. So the family is genuinely not degenerate, for two reasons: mu_rel_geo = mu_rel_helio - an Earth velocity term is a vector subtraction and so not rotationally symmetric (moving t_E), and pi_E is parallel to mu_rel, so the rotation rotates the parallax vector, which annual parallax observes.

Possible resolutions

A. Leave the gate as the resting place. Cost: the arbitrary direction persists on the gated events, and #93 is really "unresolved for literature configs". Cheapest, and defensible if the residual cost there turns out small.

B. Add mu_rel_pa as a derived symbol in mulensing/symbolic_physics.py:

mu_ra_rel  = mu_rel_mag * sin(mu_rel_pa)
mu_dec_rel = mu_rel_mag * cos(mu_rel_pa)

Then the prior's direction and the data's magnitude land on different symbols and provenance resolves them per-symbol, with the gates dropped. Note this is not a reparameterization of the sampled space -- star.pm_ra/pm_dec remain the sampled coordinates and mu_*_rel are derived -- so it changes initialization only. Bonus: it deletes the symmetric quadratic mu_rel_mag**2 = mu_ra_rel**2 + mu_dec_rel**2 that caused #92, retiring that bug class rather than sorting around it. Risks: sympy solving through sin/cos brings branch ambiguity, and adding relations perturbs the engine for every microlensing fit -- now checkable, since tests/test_seed_quality.py exists.

C. Targeted per-case handling. Derive pm from the measured pi_E when one is supplied (rather than skipping), and seed multi-source events per source rather than giving every source the same bulge mean. Smaller than B, but accumulates special cases.

D. Reparameterize the sampled coordinates -- e.g. sample mu_rel plus a mean pm instead of per-star pm, decorrelating the data-constrained direction from the prior-only directions (the same shear argument as log_q). This is the only option that could be justified on sampling grounds, which is the bar that matters; it would need ESS/Rhat evidence, measured the way the log_q work was.

Recommended next step, when this is picked up

Measure before changing anything. Run scripts/diag_rotation.py against ob161003 and ob140939 to quantify what the arbitrary direction costs on the events that are still gated. If it is small, A is the answer and this issue can close honestly. If it is large, that is the quantified case for B.

That ordering is not pedantry: twice in the work behind #99 a fix went in ahead of a measurement (seeding everything, then re-pinning a stale constant), and the measurement is what caught both.

Related gap, arguably its own issue

tests/test_seed_quality.py pins chi2/N at the seed for the four static microlensing examples, all literature configs. Nothing guards the automated MMEXOFAST path -- the one this feature is aimed at -- because it needs the Roman DC2018 data that examples/DC2018 fetches. So the path we care most about is the one measured least.

Tools now in the tree

• scripts/diag_mulens.py -- cross-machine start-value diff (versions, initvals, physical values at a fixed raw point, residual summary).
• scripts/diag_rotation.py -- the degeneracy test above; pins all four proper-motion components and aborts rather than reporting a contaminated delta.
• tests/test_seed_quality.py -- chi2/N ceilings at the seed for the four examples.

[rpoleski]

@mjmroz is working on adding parallax to MMEXOFAST. Maybe it would be best to take seed pi_E from there.

[jdeast]

I'll keep this issue open for now, but I currently think it's not really a problem at all, and generic galactic model seeds are fine.

I'm no longer convinced we need mmexofast. I've cannibalized what I need and I think i can find my own global solutions from there faster (seconds, not hours), more generally (parallax, stellar binaries, and maybe triples), and without pulling in mmexofast and sfit_minimizer as dependencies. 

I think I've basically solved this for static (no significant lens orbital motion) 2L1S systems. I'm waiting in the results from a sweep of all 30 such systems in the 2018 data challenge, but smaller subset (5/5) worked to recover observable and physical parameters from the truth table.