commit 857b086fcddd4bb1476bc5a5892b27e3bb4ccead from: ale date: Fri Sep 11 03:25:39 2026 UTC Add digest2 NEO scoring and cross-night tracklet linking Investigated 4 candidate technologies for improving the pipeline; GPU reprojection and local plate-solving were rejected with real reasons (no GPU path in Reproject.jl/wcslib; multi-GB index files for marginal gain), documented in Design refinements. Built the other two: - digest2_score wraps the MPC's own external digest2 classifier (verified by actually cloning, building, and running it). Real verification against ZTF field 451 surfaced a genuinely useful, non-obvious finding: it correctly scored the field's 2 known Main Belt objects low on neo_score, while 131 of 133 tracklets scored 100 — bogus links of stationary stars, an artifact of that demo's own loose match_radius, not something digest2 got wrong. Written up in full on its own wiki page. - link_across_nights extends link_candidates' within-night linear linking across multiple observing nights (fit-and-extrapolate on sky coordinates, union-find grouping). Validated synthetically; a real multi-night same-object case was searched for (checked the existing 5-field IASC/PS1 dataset, then queried SkyBoT directly) but not found within reasonable effort — stated plainly as an open gap rather than skipped past. 188/188 tests pass (digest2_score's live test ran for real against a locally-built binary, not just skipped). commit - 1acf0bc70592e4401c60ae2d4b010f082fea8205 commit + 857b086fcddd4bb1476bc5a5892b27e3bb4ccead blob - 1d0f8ba7d8154d76cfa9c4122819c596a10a35ba blob + 7d069a471f88881980e38a7dded92bf7d8193b00 --- .env.example +++ .env.example @@ -5,3 +5,9 @@ # at https://nova.astrometry.net/. Enables the live plate_solve round-trip # test; without it, that one test is skipped instead of run. ASTROMETRY_API_KEY= + +# Path to a locally-built digest2 executable (see docs/src/api/digest2.md +# for build steps). Enables the live digest2_score test; without it, +# that one test is skipped instead of run. Defaults to "digest2" (PATH +# lookup) if unset. +DIGEST2_PATH= blob - ac796e08ccc95aa709ee8492096cdd77772d876a blob + 57bbddc167d5d13fdccdf1042a1464b7a204e0f5 --- docs/make.jl +++ docs/make.jl @@ -64,6 +64,8 @@ makedocs(; "Pipeline" => "api/pipeline.md", "Rotation Period" => "api/rotation.md", "MPC / ADES Export" => "api/mpc-export.md", + "MPC digest2 Scoring" => "api/digest2.md", + "Cross-Night Linking" => "api/cross-night-linking.md", ], ], ) blob - /dev/null blob + 5063549e404246fb1f4b853dbee10e4779b4c7ac (mode 644) --- /dev/null +++ docs/src/api/cross-night-linking.md @@ -0,0 +1,83 @@ +# Cross-Night Linking + +[`link_candidates`](@ref) only links detections *within* one observing +night (a flat sequence of frames sharing one linear-motion model over +minutes to hours). [`link_across_nights`](@ref) solves the same kind of +problem one level up: grouping different nights' own tracklets that are +consistent with one real object, over gaps of days where real orbital +curvature — not just the instantaneous rate — starts to matter. This was +entirely greenfield: nothing in this pipeline had any night/date/session +concept before this (confirmed by a direct search of `src/` — `timestamps` +was, and still is for a single night, one flat `Vector{Float64}` of +Julian Dates in lockstep with `fits_paths`). + +## Example + +```julia +night1 = search_field(night1_paths; timestamp_key="OBSMJD").movers +night2 = search_field(night2_paths; timestamp_key="OBSMJD").movers +night3 = search_field(night3_paths; timestamp_key="OBSMJD").movers + +groups = link_across_nights([night1, night2, night3]) +``` + +returns a `Vector` of groups, each a `Vector` of `(night, id)` pairs — +tracklets, from different nights, judged to be the same real object. +Only groups spanning at least `min_nights` (default `2`) distinct nights +are returned. + +## Design + +Each night's tracklets are independently fit to a linear rate in a +local tangent-plane projection (reusing `_linfit` — the exact +same fitting routine [`link_candidates`](@ref)'s own within-night refit +already uses, not reimplemented), then extrapolated forward to every +other night's tracklets' own mean epoch. A pair is accepted within +`match_radius_arcsec` of the extrapolated position, optionally loosened +by `max_accel_arcsec_per_day2 * Δt_days^2` to allow for real curvature +over longer gaps (off, `Inf`, by default — mirroring +[`link_candidates`](@ref)'s own `max_speed::Real=Inf` default: a single +generous radius does the real work unless the caller has independent +reason to tighten it). Accepted pairs are grouped transitively via +union-find, since one tracklet can plausibly pair with several other +nights' tracklets that must all collapse into one group. + +## An honest gap: not yet validated against a real multi-night case + +Every other real-data claim in this project's docs is backed by an +actual run against real data. This one isn't yet, and that's stated +plainly rather than skipped past. + +The existing real IASC/PS1 dataset was checked first (5 fields, 9 known +objects — see [Validating against real IASC (Pan-STARRS1) campaign data](../iasc-campaign-validation.md)): none of +those 9 objects repeats across fields, so there was no ready-made real +cross-night case sitting in data already on hand. A further, bounded +real check — querying SkyBoT for ZTF field 451's own sky position at +its original epoch and again 3 and 7 days later — confirmed *why* this +is genuinely hard to find by chance: every object SkyBoT reports within +that field changes completely from one query to the next. Even a +"slow" Main Belt object here (2002 UY45, ~852"/day) crosses a ZTF +quadrant's own ~35' width in a few days — real multi-night, same-tile +recovery needs an object caught unusually close to its apparent +stationary point, not just "a slow-ish asteroid," and finding one that +also happens to fall within ZTF's actual public observing history for +one specific tile is a real needle-in-a-haystack search, not something +a couple of cone searches turns up. + +So: [`link_across_nights`](@ref)'s correctness is currently established +by synthetic exact-recovery tests only (construct a real linear-motion +object across 3 fabricated nights, plus unrelated same-night "noise" +tracklets, and confirm the object's tracklets group correctly while +noise never does — mirroring [`link_candidates`](@ref)'s own synthetic +test style exactly). `match_radius_arcsec`'s default (5") and +`max_accel_arcsec_per_day2=Inf` are starting points to tune against a +real campaign's own cadence and target population, not values +calibrated against this project's own measured data the way e.g. +`variability_chi2`'s `systematic_error_fraction` was. A real multi-night +validation remains open for whenever a suitable real case turns up (a +live campaign's own multi-night data would be the natural source). + +```@autodocs +Modules = [AsteroidPipeline] +Pages = ["linking_multinight.jl"] +``` blob - /dev/null blob + a7f1796b284da03164bc58d8c2bd6493eea7a343 (mode 644) --- /dev/null +++ docs/src/api/digest2.md @@ -0,0 +1,86 @@ +# MPC digest2 Scoring + +Scores each tracklet in an [`astrometric_calibrate`](@ref)-shaped table +for its likelihood of belonging to several real solar-system orbit +classes (NEO, Main Belt, Jupiter Trojan, ...), via `digest2` — the +Minor Planet Center's own real, already-validated short-arc classifier +(statistical ranging against a real population model), not anything +trained or guessed at here. Explored alongside three other candidate +technologies (see [Design refinements](@ref) for why GPU-accelerated +reprojection and local plate-solving were investigated and *not* +recommended); `digest2` and [Cross-Night Linking](@ref) were the two +worth building. + +## Setup (one-time, outside this package) + +`digest2` is a separate C program, not a Julia dependency: + +```sh +git clone https://github.com/Smithsonian/digest2.git +cd digest2/digest2 +make +``` + +The `MPC.config`, `digest2.model.csv`, and `digest2.obscodes` files this +needs at runtime already ship in that same `digest2/` directory — +nothing to write or copy yourself. Either put the resulting `digest2` +executable on `PATH`, or pass its path directly: + +```julia +digest2_score(candidates, "I41"; digest2_path="/path/to/digest2/digest2/digest2") +``` + +## Example + +```julia +scores = digest2_score(candidates, "I41") +``` + +produces one row per tracklet: + +| id | rms | int_score | neo_score | n22_score | n18_score | +|--:|--:|--:|--:|--:|--:| +| 1 | 0.02 | 100.0 | 100.0 | 23.0 | 0.0 | +| 2 | 0.0 | 100.0 | 98.0 | 26.0 | 1.0 | + +## A real finding: this is an orbit-class classifier, not a real/bogus one + +Run for real against `real_data_demo.jl`'s ZTF field 451 baseline (133 +tracklets from real, already-downloaded ZTF data; setup: `examples/fetch_data.sh` +then `examples/real_data_demo.jl`'s baseline stage), the result was the +*opposite* of a naive expectation: + +| | `neo_score` | +|:--|--:| +| The 2 real, SkyBoT-confirmed known objects (2002 UY45, 1997 KO3 — both Main Belt) | 5, 7 | +| The other 131 tracklets | 100 (130 of 131) | + +That's `digest2` working correctly, not a bug: it correctly recognized +the 2 real objects as *not* NEOs (they're Main Belt — their score is +concentrated in `digest2`'s `MB1`/`MB2` classes, which this function +doesn't parse out). The other 131 were never real moving objects at +all — inspecting one directly, its 5 detections jitter by under 1" +across the full 6.25 h baseline (`~3"/day` implied rate): a real, +stationary star, re-detected each frame and linked into a bogus +tracklet only because [`link_candidates`](@ref)'s `match_radius` in +that demo (10", looser than ZTF's real sub-arcsec centroiding +precision — the same lesson [Validating against real IASC (Pan-STARRS1) campaign data](../iasc-campaign-validation.md) +already learned once from `PERROR`, resurfacing here) was loose enough +to accept a star's own centroid jitter as if it were motion. `digest2` +has no way to tell a stationary star's jitter from a real, very distant +object moving at a dynamically-consistent near-zero rate — it scored +the hypothesis it was given correctly; the hypothesis itself was never +a real tracklet. + +**Practical takeaway**: don't sort by `neo_score` on raw +[`link_candidates`](@ref) output as a "most interesting first" filter. +Tighten `match_radius` to the survey's own real astrometric precision +first, or otherwise filter for genuinely consistent motion, before +scoring — `digest2_score` classifies an *already-plausible* tracklet's +dynamical class well; it is not a substitute for that upstream quality +control. + +```@autodocs +Modules = [AsteroidPipeline] +Pages = ["digest2.jl"] +``` blob - c29d9e2f32b246d5f916e654ca09cd2ee294589a blob + bdee5ac012ec934eeb2efae0227a6cc2fe19a7c0 --- docs/src/design-refinements.md +++ docs/src/design-refinements.md @@ -279,3 +279,43 @@ this profiling pass found reason to change. No second `build_reference`-sized win turned up. The pipeline's real remaining cost, at typical real field densities, is dominated by what was already found and fixed. + +## Four candidate technologies investigated; two built, two rejected with real reasons + +Asked generally what other technology could improve this pipeline +(separately from whether a boosted-decision-tree real/bogus classifier +was worth adding — rejected on its own, before this: only ~12 real +confirmed positives exist across every real run this project has done, +nowhere near enough to train or validate one). Four concrete candidates, +each checked against this project's own real environment/code rather +than decided on general reputation: + +- **GPU-accelerated reprojection** (`build_reference`'s real bottleneck, + already solved 3.21x via multiprocessing) — rejected. `Reproject.jl`'s + own source has zero GPU code path, and the actual per-pixel cost is + `wcslib`'s C calls, which aren't portable to a GPU kernel without + reimplementing WCS math from scratch — real risk in a domain that has + already produced two real bugs here (the transposed-aperture bug, the + cross-process `WCSTransform` segfault this same multiprocessing work + hit), for a target already solved. +- **Local plate-solving** (replacing `plate_solve`'s live + nova.astrometry.net dependency) — rejected. `solve-field` isn't in + nixpkgs, and needs multiple GB of scale-specific index files even + once built, to replace something that already works. +- **`digest2` NEO/orbit-class scoring** — built (`digest2_score`, [MPC + digest2 Scoring](@ref)). Real verification against real ZTF field 451 + candidates surfaced a genuinely useful, non-obvious finding: `digest2` + correctly scored the field's 2 real known objects (Main Belt) low on + `neo_score`, while 131 of the other 133 tracklets — bogus links of + ordinary stationary stars, an artifact of that demo's own + looser-than-ZTF's-real-precision `match_radius` — scored `neo_score=100`. + Not a bug in `digest2_score`; a real lesson about what it can and + can't tell apart, written up in full on its own page. +- **Cross-night tracklet linking** — built (`link_across_nights`, + [Cross-Night Linking](@ref)). Genuinely greenfield (no night/session + concept existed anywhere in `src/` before this). Validated + synthetically only so far — a real check (querying SkyBoT for field + 451's own sky position 3 and 7 days out) confirmed *why* a real + multi-night same-object case is hard to find by chance: every object + in the field changes completely within days, since even a "slow" + Main Belt object outpaces a ZTF quadrant's own width in under a week. blob - 17d1bcf1ff5643cd933349eb71333d684b8c8853 blob + 467e48622b674c84e89c1c1ee9267d2f7f3a873f --- docs/src/index.md +++ docs/src/index.md @@ -53,6 +53,13 @@ does the same for the legacy fixed-width 80-column for required by some programs (e.g. IASC, as of 2026) even though the MPC's own submissions now prefer ADES. +`digest2_score` scores tracklets against real solar-system orbit +classes via the Minor Planet Center's own external `digest2` classifier +— see [MPC digest2 Scoring](@ref), including a real finding about what +it can and can't tell apart. `link_across_nights` extends +`link_candidates`'s within-night linking across multiple observing +nights — see [Cross-Night Linking](@ref). + ## Status !!! note "Early development, but validated end to end against real data" @@ -69,6 +76,8 @@ own submissions now prefer ADES. | `fit_moffat_psf` (PSF analytic fallback) | Synthetic · PSF width calibrated on real ZTF data | | `plate_solve` | Live nova.astrometry.net service | | `crossmatch_catalog` (`:skybot`/`:vsx`/`:simbad`) | Live services, real positive controls | +| `digest2_score` | Real, locally-built `digest2` binary · real ZTF field 451 candidates (see [MPC digest2 Scoring](@ref) for a real, non-obvious finding) | +| `link_across_nights` | Synthetic exact-recovery only — no real multi-night same-object case found yet, see [Cross-Night Linking](@ref) | Not yet run on a live IASC search campaign (as opposed to practice data) — see [Using real IASC campaign data](@ref) for what that would involve. @@ -262,6 +271,24 @@ real campaign, or any other survey's data: a candidate as needing independent confirmation (a catalog match or a recovered period), not as self-evidently real. +!!! note "GPU reprojection and local plate-solving: investigated, not recommended now" + Explored alongside `digest2`/cross-night linking as candidate + technologies. **GPU-accelerated reprojection**: `Reproject.jl` (used + by `build_reference`) has no GPU code path at all (checked its + source directly), and the real per-frame cost is `wcslib`'s own C + calls — not portable to a GPU kernel without reimplementing WCS + pixel math from scratch, real correctness risk in a domain that has + already produced two real bugs this project (the transposed-aperture + bug, the cross-process `WCSTransform` segfault), for a target + already solved by a real, measured 3.21x via multiprocessing (see + [Design refinements](@ref)). **Local plate-solving** (`solve-field`, + replacing the live nova.astrometry.net API `plate_solve` already + uses): not packaged in nixpkgs, and needs multiple GB of scale-specific + index files even once built — real setup cost to replace something + that already works. Recorded here so this isn't silently forgotten + or re-litigated from scratch without new information changing the + calculus. + ## Follow-up workflows ### Rotation period recovery blob - a3264f088ab0d4d2124bd3df972cdcd3cef28e6c blob + 9acc10405471cacc11bc605f28edb8f25170aeb7 --- src/AsteroidPipeline.jl +++ src/AsteroidPipeline.jl @@ -18,6 +18,7 @@ using Distributed include("detection.jl") include("linking.jl") +include("linking_multinight.jl") include("variables.jl") include("astrometry.jl") include("platesolve.jl") @@ -28,11 +29,12 @@ include("zogy.jl") include("pipeline.jl") include("rotation.jl") include("mpc_export.jl") +include("digest2.jl") export detect_sources, link_candidates, load_wcs, pix_to_sky, astrometric_calibrate, crossmatch_catalog, run_pipeline, build_reference, load_frame, estimate_psf, fit_moffat_psf, zogy_subtract, light_curve, recover_rotation_period, plate_solve, search_field, find_variable_sources, variability_chi2, photometric_scale, - ades_psv, mpc80_report, julian_date_to_iso8601 + ades_psv, mpc80_report, julian_date_to_iso8601, digest2_score, link_across_nights end # module AsteroidPipeline blob - /dev/null blob + cef600f602d7e26238aed613090f85d9a009dd5f (mode 644) --- /dev/null +++ src/digest2.jl @@ -0,0 +1,117 @@ +""" + digest2_score(candidates, station::AbstractString; + digest2_path::AbstractString="digest2", + config_dir::Union{Nothing,AbstractString}=nothing, + trksub_prefix::AbstractString="") -> Table + +Score each tracklet in `candidates` (an [`astrometric_calibrate`](@ref) +table, same shape [`mpc80_report`](@ref) takes) for its likelihood of +belonging to several real orbit classes, via the Minor Planet Center's +own `digest2` — a real, external, already-validated classifier (short-arc +statistical ranging against a real solar-system population model), not +anything trained or guessed at here. + +`digest2` is a separate C program (source: +`https://github.com/Smithsonian/digest2`), not a Julia dependency — +build it yourself (`make` in its `digest2/` directory; see +[MPC digest2 Scoring](@ref) for exact steps) and either put the +resulting `digest2` executable on `PATH` or pass its path via +`digest2_path`. `config_dir` defaults to `dirname(digest2_path)` (its own +documented "same directory as the executable" layout) and must contain +`digest2.model.csv`, `digest2.obscodes`, and an `MPC.config` — all three +ship with the `digest2` repository itself and are used unmodified here +(this pipeline supplies none of its own: `digest2`'s upstream `MPC.config` +already sets a real, tuned per-observatory error, including ZTF's own +`I41`, and rewriting it would just be guessing at numbers `digest2`'s own +maintainers already calibrated). + +Reuses [`mpc80_report`](@ref) verbatim to build `digest2`'s input (`digest2` +requires MPC 80-column input grouped by consecutive same-designation +lines, sorted by designation then time — `candidates` is sorted by +`(id, epoch)` locally first, since [`astrometric_calibrate`](@ref) only +guarantees grouping by `id`, not epoch order across the whole table, and +a caller's own row order should not silently change `digest2`'s result). +Piped via `stdin`/`stdout` through a `Cmd` built from an argument vector +(never a shell string), so `digest2_path`/`config_dir` cannot be +mis-parsed as shell syntax. + +Returns one row per tracklet (`id`, matching `candidates`' own `id`) with +`digest2`'s `RMS` (its own linear-motion fit residual, arcsec — a +free byproduct neither [`link_candidates`](@ref) nor +[`astrometric_calibrate`](@ref) compute) and four raw 0-100 +pseudo-probability scores: `int_score` (any orbit class of general MPC +interest), `neo_score` (near-Earth object), `n22_score`/`n18_score` (NEO +with absolute magnitude ≤ 22/≤ 18 respectively — i.e. how large the +object would have to be). These are `digest2`'s own default output +columns, present regardless of `MPC.config`'s requested class list (real +behavior, confirmed by actually building and running `digest2`, not +assumed from its docs) — this function does not parse the trailing +"Other Possibilities" free-text column (main-belt/Trojan/comet class +scores for whatever doesn't make the main four). + +**Real finding, not a design assumption — `digest2`'s scores are an +orbit-class classifier, not a real/bogus one, and this matters in +practice**: run for real against `real_data_demo.jl`'s ZTF field 451 +baseline tracklets (see [MPC digest2 Scoring](@ref) for the full +numbers), the two real, SkyBoT-confirmed known objects scored `neo_score` +5 and 7 (correctly low — they're Main Belt, not NEOs), while 131 of the +other 133 tracklets scored `neo_score=100`. Those 131 were not missed +NEOs: inspecting one directly, its 5 points jitter by under 1" across +the full 6.25 h baseline (`~3"/day` implied rate — this is a real, +stationary star, re-detected each frame and "linked" into a bogus +tracklet only because [`link_candidates`](@ref)'s `match_radius` here +(10", well above ZTF's real sub-arcsec centroiding precision) is loose +enough to also accept a star's own centroid jitter as if it were +motion). `digest2` has no way to know that near-zero, jittery apparent +motion is a stationary star rather than a real object at a range where +that motion is dynamically self-consistent — it scored the *hypothesis* +correctly, on input that was never a real tracklet to begin with. + +**The actionable lesson**: don't sort by `neo_score` on raw +[`link_candidates`](@ref) output as a "most interesting first" filter — +tighten `match_radius` to the survey's own real astrometric precision +first (`iasc-campaign-validation.md` already established this exact +principle from PS1's `PERROR`), or otherwise vet tracklets for genuine, +consistent motion, before scoring; `digest2_score` is a real, valuable +tool for classifying the dynamical class of an *already-plausible* +tracklet, not a substitute for that upstream quality control. + +Throws `ArgumentError` if `digest2_path` resolves to nothing runnable +(`Sys.which`-style — an explicit argument, not an environment variable +this function reads itself, matching [`plate_solve`](@ref)'s own +`api_key` convention) or if `station` isn't a 3-character MPC +observatory code. +""" +function digest2_score(candidates, station::AbstractString; + digest2_path::AbstractString="digest2", + config_dir::Union{Nothing,AbstractString}=nothing, + trksub_prefix::AbstractString="") + length(station) == 3 || throw(ArgumentError("station must be a 3-character MPC observatory code")) + + resolved = isfile(digest2_path) ? digest2_path : Sys.which(digest2_path) + resolved === nothing && throw(ArgumentError( + "digest2 executable not found at or on PATH as \"$digest2_path\" — build it " * + "from https://github.com/Smithsonian/digest2 and pass its path via digest2_path")) + dir = config_dir === nothing ? dirname(abspath(resolved)) : config_dir + + sorted = sort(candidates; by=row -> (row.id, row.epoch)) + input = mpc80_report(sorted, station; trksub_prefix=trksub_prefix) + + out = IOBuffer() + run(pipeline(`$resolved -p $dir -`, stdin=IOBuffer(input), stdout=out)) + lines = split(chomp(String(take!(out))), "\n") + + ids, rmss, ints, neos, n22s, n18s = Int[], Float64[], Float64[], Float64[], Float64[], Float64[] + for line in lines[2:end] # skip the "Desig. RMS Int NEO N22 N18 ..." header + tokens = split(line) + designation = tokens[1][(length(trksub_prefix)+1):end] + push!(ids, parse(Int, designation; base=36)) + push!(rmss, parse(Float64, tokens[2])) + push!(ints, parse(Float64, tokens[3])) + push!(neos, parse(Float64, tokens[4])) + push!(n22s, parse(Float64, tokens[5])) + push!(n18s, parse(Float64, tokens[6])) + end + + return Table(id=ids, rms=rmss, int_score=ints, neo_score=neos, n22_score=n22s, n18_score=n18s) +end blob - /dev/null blob + a74ba256ab773f60d0a12fb02fd56971ec8b8614 (mode 644) --- /dev/null +++ src/linking_multinight.jl @@ -0,0 +1,109 @@ +""" + link_across_nights(night_movers::AbstractVector; + match_radius_arcsec::Real=5.0, + max_accel_arcsec_per_day2::Real=Inf, + min_nights::Integer=2) + +Group tracklets from *different* observing nights that are consistent +with one real moving object, the same problem [`link_candidates`](@ref) +solves within a single night, one level up: on sky coordinates rather +than pixels (different nights can have entirely different pointings/WCS +solutions), and over gaps of days rather than minutes, where the target's +own real orbital curvature — not just its instantaneous rate — starts to +matter. + +`night_movers[k]` is night `k`'s own [`astrometric_calibrate`](@ref) (or +[`search_field`](@ref)) output table (`id`, `frame`, `x`, `y`, `ra`, +`dec`, `epoch`; `x`/`y` are unused here). Each night's tracklets (per +`id`) are independently fit to a linear rate in a local tangent-plane +projection around their own mean position — via [`_linfit`](@ref), +reused as-is from [`link_candidates`](@ref)'s own within-night fit, not +reimplemented — then that rate is used to extrapolate forward to every +other night's tracklets' own mean epoch. A cross-night pair is accepted +whenever the extrapolated position lands within `match_radius_arcsec` of +the other tracklet's actual mean position, loosened (if +`max_accel_arcsec_per_day2` is finite) by an added +`max_accel_arcsec_per_day2 * Δt_days^2` allowance for real curvature +over the longer gap. Accepted pairs are grouped transitively (one +tracklet can plausibly pair with several other nights' tracklets, which +must then all collapse into a single group) via union-find; only groups +spanning at least `min_nights` distinct nights are kept. Tracklets with +fewer than 2 points (no rate to fit) are skipped — they carry no motion +information to extrapolate from and so cannot be matched here. + +`max_accel_arcsec_per_day2` defaults to `Inf` (off), mirroring +[`link_candidates`](@ref)'s own `max_speed::Real=Inf` default exactly: +a single generous `match_radius_arcsec` does the real matching work by +default. **Neither default is derived from this project's own measured +data** — unlike e.g. `variability_chi2`'s `systematic_error_fraction` +(calibrated against three real confirmed variables), no real multi-night +same-object case was available to calibrate against when this was +written; treat `match_radius_arcsec=5.0` as a starting point to tune +against your own campaign's real cadence and target population, not a +validated constant. + +Returns a `Vector` of groups, each a `Vector` of `(night, id)` +`NamedTuple`s — the tracklets, across nights, judged to be the same +object — mirroring [`link_candidates`](@ref)'s own "vector of vectors" +return shape rather than a flat table. +""" +function link_across_nights(night_movers::AbstractVector; + match_radius_arcsec::Real=5.0, + max_accel_arcsec_per_day2::Real=Inf, + min_nights::Integer=2) + Node = NamedTuple{(:night, :id),Tuple{Int,Int}} + nodes = Node[] + fits = NamedTuple[] # (ra0, dec0, epoch0, rate_ra, rate_dec) per node, same order as nodes + + for (night, candidates) in enumerate(night_movers) + for id in unique(candidates.id) + rows = collect(filter(r -> r.id == id, candidates)) + length(rows) < 2 && continue + + ra0, dec0 = mean(r.ra for r in rows), mean(r.dec for r in rows) + epoch0 = mean(r.epoch for r in rows) + cos_dec0 = cosd(dec0) + t = [r.epoch - epoch0 for r in rows] + dra = [(r.ra - ra0) * cos_dec0 * 3600 for r in rows] + ddec = [(r.dec - dec0) * 3600 for r in rows] + ra_int, rate_ra = _linfit(t, dra) + dec_int, rate_dec = _linfit(t, ddec) + + push!(nodes, (night=night, id=id)) + push!(fits, (ra0=ra0 + ra_int / cos_dec0 / 3600, dec0=dec0 + dec_int / 3600, + epoch0=epoch0, rate_ra=rate_ra, rate_dec=rate_dec, cos_dec0=cos_dec0)) + end + end + + n = length(nodes) + parent = collect(1:n) + find(i) = (while parent[i] != i; i = parent[i]; end; i) + function union!(i, j) + ri, rj = find(i), find(j) + ri != rj && (parent[ri] = rj) + end + + for i in 1:n, j in 1:n + i == j && continue + nodes[i].night == nodes[j].night && continue + a, b = fits[i], fits[j] + dt = b.epoch0 - a.epoch0 + dt <= 0 && continue # only extrapolate forward in time; the (j,i) pair covers the reverse + + pred_ra = a.ra0 + (a.rate_ra * dt) / a.cos_dec0 / 3600 + pred_dec = a.dec0 + a.rate_dec * dt / 3600 + sep_arcsec = hypot((pred_ra - b.ra0) * a.cos_dec0, pred_dec - b.dec0) * 3600 + + tolerance = match_radius_arcsec + isfinite(max_accel_arcsec_per_day2) && (tolerance += max_accel_arcsec_per_day2 * dt^2) + + sep_arcsec <= tolerance && union!(i, j) + end + + groups = Dict{Int,Vector{Node}}() + for i in 1:n + push!(get!(groups, find(i), Node[]), nodes[i]) + end + + return [g for g in values(groups) if length(unique(node.night for node in g)) >= min_nights] +end blob - 951f8032889d36093149f46cd20a5de048b3b2a9 blob + 5a9b7fa0f2e6d3ccdf4ef12e78b0d21d188ed1c0 --- test/runtests.jl +++ test/runtests.jl @@ -1140,4 +1140,96 @@ end @test_throws ArgumentError mpc80_report(candidates, "I41"; trksub_prefix="TOOLONG") end + @testset "digest2_score" begin + # digest2 is a separate, external C program (not a Julia + # dependency) — never built as part of Pkg.test(); skip cleanly + # when it isn't installed, same pattern as plate_solve's + # ASTROMETRY_API_KEY skip below. + digest2_path = get(ENV, "DIGEST2_PATH", "digest2") + resolved = isfile(digest2_path) ? digest2_path : Sys.which(digest2_path) + + candidates = Table(id=[1, 1, 1], frame=[1, 2, 3], x=[1.0, 2.0, 3.0], y=[1.0, 2.0, 3.0], + ra=[150.0, 150.05, 150.10], dec=[20.0, 20.02, 20.04], + epoch=[2451545.0, 2451545.04, 2451545.08]) + + @test_throws ArgumentError digest2_score(candidates, "XX") # not 3 characters + @test_throws ArgumentError digest2_score(candidates, "I41"; digest2_path="not-a-real-digest2-binary") + + if resolved === nothing + @test_skip "digest2 not installed" + else + scores = digest2_score(candidates, "I41"; digest2_path=resolved) + @test length(scores) == 1 + @test scores.id[1] == 1 + @test 0.0 <= scores.neo_score[1] <= 100.0 # a real score, not this test asserting a specific value + end + end + + @testset "link_across_nights" begin + T0 = 2451545.0 + RATE_RA, RATE_DEC = 33.83, 18.0 # arcsec/day, in the tangent plane around (150,20) + + # a real object: exactly linear motion (in the tangent plane), + # split across 3 nights days apart. Night 2 and 3 continue the + # same global rate, so an independent per-night fit recovers it + # exactly and the cross-night extrapolation has zero residual. + function object_row(t) + dt = t - T0 + dec = 20.0 + RATE_DEC * dt / 3600 + ra = 150.0 + RATE_RA * dt / (cosd(20.0) * 3600) + return (ra=ra, dec=dec, epoch=t) + end + night1_obj = [object_row(T0 - 0.05), object_row(T0 + 0.05)] + night2_obj = [object_row(T0 + 4.95), object_row(T0 + 5.05)] + night3_obj = [object_row(T0 + 9.95), object_row(T0 + 10.05)] + + # "noise" tracklets: real, self-consistent motion within their own + # night (so link_across_nights can fit a rate at all), but at an + # unrelated sky position/rate that must never match the real + # object's extrapolated position across nights. + noise(t0, ra0, dec0) = [(ra=ra0, dec=dec0, epoch=t0), (ra=ra0 + 0.01, dec=dec0 - 0.01, epoch=t0 + 0.1)] + night1_noise = noise(T0, 200.0, -10.0) + night2_noise = noise(T0 + 5.0, 44.0, 60.0) + night3_noise = noise(T0 + 10.0, 300.0, 0.0) + + to_table(obj, noise) = Table( + id=vcat(fill(1, length(obj)), fill(2, length(noise))), + frame=vcat(1:length(obj), 1:length(noise)), + x=zeros(length(obj) + length(noise)), y=zeros(length(obj) + length(noise)), + ra=vcat([r.ra for r in obj], [r.ra for r in noise]), + dec=vcat([r.dec for r in obj], [r.dec for r in noise]), + epoch=vcat([r.epoch for r in obj], [r.epoch for r in noise])) + + night_movers = [to_table(night1_obj, night1_noise), to_table(night2_obj, night2_noise), + to_table(night3_obj, night3_noise)] + + groups = link_across_nights(night_movers) + @test length(groups) == 1 # only the real object spans >= min_nights=2 + real_group = groups[1] + @test length(real_group) == 3 # one tracklet per night + @test Set(g.night for g in real_group) == Set([1, 2, 3]) + @test all(g.id == 1 for g in real_group) # never picks up a noise id=2 tracklet + + # min_nights: raising it above what the real object spans drops it too + @test isempty(link_across_nights(night_movers; min_nights=4)) + + # a real 3" offset from the exact extrapolation is still within + # the default 5" tolerance, but not a tightened 1" one + offset_row(t) = merge(object_row(t), (dec=object_row(t).dec + 3.0 / 3600,)) + night2_offset = [offset_row(T0 + 4.95), offset_row(T0 + 5.05)] + offset_movers = [night_movers[1], to_table(night2_offset, night2_noise), night_movers[3]] + + offset_groups = link_across_nights(offset_movers) + @test length(offset_groups) == 1 + @test length(offset_groups[1]) == 3 # still spans all 3 nights within the default tolerance + + # tightened past the 3" offset: night 2's tracklet drops out, but + # nights 1 and 3 (never offset, and 10 days apart is still an exact + # extrapolation under this exactly-linear synthetic model) still + # correctly link on their own + tight_groups = link_across_nights(offset_movers; match_radius_arcsec=1.0) + @test length(tight_groups) == 1 + @test Set(g.night for g in tight_groups[1]) == Set([1, 3]) + end + end