commit - 1acf0bc70592e4401c60ae2d4b010f082fea8205
commit + 857b086fcddd4bb1476bc5a5892b27e3bb4ccead
blob - 1d0f8ba7d8154d76cfa9c4122819c596a10a35ba
blob + 7d069a471f88881980e38a7dded92bf7d8193b00
--- .env.example
+++ .env.example
# 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
"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
+# 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
+# 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
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
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"
| `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.
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
include("detection.jl")
include("linking.jl")
+include("linking_multinight.jl")
include("variables.jl")
include("astrometry.jl")
include("platesolve.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
+"""
+ 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
+"""
+ 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
@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
