Smithsonian / Smithsonian/layup

The weight matrix ignores `Observation::inverse_covariance`, so an astrometric uncertainty cannot be oriented

Open
#583 1 comment 0 reactions 0 assignees View on GitHub

Nobody has claimed this yet.

Dominant language
Python
Stars
10
Forks
2
Avg merge
1d 43m
Merged PRs (30d)
58

Description

Observation carries an optional inverse covariance (detection.cpp:150, exposed to Python at detection.cpp:413), but the fitter never reads it. get_weight_matrix builds W from ra_unc and dec_unc alone:

orbit_fit.cpp:1034   double x_unc = *detections[i].ra_unc;
                     tripletList.push_back(T(r0,   r0,   1.0 / (x_unc * x_unc)));
                     double y_unc = *detections[i].dec_unc;
                     tripletList.push_back(T(r0+1, r0+1, 1.0 / (y_unc * y_unc)));

Two diagonal entries and no off-diagonal, so the uncertainty of an astrometric observation is always an ellipse aligned with the right ascension and declination axes. The two axes may differ, so it is not isotropic, but it cannot be rotated.

That is the wrong shape for one real effect. A timing error puts the positional error along the direction of the object's angular motion described in Farnocchia et al. (2015), so the correct uncertainty is an ellipse elongated along the sky-motion direction, which in general lies along neither axis. Representing it needs the off-diagonal term, which is what inverse_covariance is for.

The effect is not small for fast movers. Scoring observations against JPL orbits on a set of near-Earth objects, with sigma replaced by sqrt(sigma^2 + (rate * sigma_t)^2) and sigma_t = 2.3 s, removes 68 per cent of the apparent outliers on objects moving faster than 2 arcsec/min (25 to 8) while moving none of the slow control (4 to 4). The 2.3 s is measured rather than fitted: residual against sky rate has a slope of 0.039 arcsec per arcsec/min, and it is consistent with the 0.864 s date quantum of the obs80 format.

That measurement had to inflate both axes, because scalar sigma is all the interface accepts, which over-inflates the cross-track direction and makes the result an upper bound rather than a demonstration. Honouring inverse_covariance when it is present would make the correct test possible, and is a contained change: a 2x2 block in place of the two diagonal entries, falling back to the present behaviour when the field is absent.

Worth deciding at the same time whether the rate-aware sigma should be computed inside Layup from the fitted sky motion, or supplied by the caller.

Drafted with Claude Opus 5.

Contributor guide

Open the contributing guide

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

Research direction

Start with get_weight_matrix in orbit_fit.cpp:1034 and the Observation definitions and Python exposure in detection.cpp:150 and :413. Trace how inverse_covariance is represented and decide how the unresolved rate-aware sigma belongs in Layup or the caller. Done means using a 2x2 covariance block when present while preserving the current diagonal fallback when it is absent.

Written by the indexing model from the issue text.

Assessment

Tech stack
cpp, python
Domain
backend, data
Issue type
Feature
Difficulty
4/5
Estimated time
3-5 days
Activity status
Active
Clarity
Mostly clear
Newbie friendliness
52/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.