commit dd2afcad8fb31f92714d78933a8e2aa88a22de0a from: ale date: Sun Aug 16 16:12:39 2026 UTC Fix detect_sources's transposed aperture bug; migrate cross-match off X-Match; recalibrate variable-source detection against real data detect_sources built CircularAperture(row.x, row.y, r) directly from PeakMesh's output, but Photometry.jl's CircularAperture indexes the array in the opposite (x=1st dim, y=2nd) sense internally from PeakMesh's own (x=2nd dim, y=1st) convention -- every measured flux was silently wrong (~140x too small in a clean test case) except when a source's row and column happened to coincide. Positions were never affected, so every prior position-based result (tracklet linking, WCS calibration, the real ZTF baseline/ZOGY tracklet counts) stands unchanged; only flux/flux_err were wrong, since nothing depended on them quantitatively until find_variable_sources. Found and fixed by comparing measured flux against the analytic enclosed-energy integral for an isolated synthetic source. That fix invalidated the real-data numbers find_variable_sources's constants were originally calibrated against (a false ~29% median flux_err/flux, a false ~2-of-120 S/N pass rate, a false selection-bias explanation for a photometric-scale mismatch that's actually ordinary seeing-dependent aperture correction) -- re-measured everything against the same real ZTF frames post-fix and rewrote photometric_scale's and find_variable_sources's docstrings and chi2_threshold default (3 -> 50) against the corrected measurements, including the honestly-reported residual false-positive rate that remains. Also: crossmatch_catalog's :vsx/:simbad backend moved from the CDS X-Match service (extended, total outages) to direct SIMBAD/VizieR TAP queries; detect_sources gained an optional gain parameter for Poisson-aware flux_err; photometric_scale added as a raw-path quality signal in run_pipeline, honestly documented as unvalidated against a known-real anomaly (unlike quality_max_std). Full test suite green (118 pass, 1 broken -- the live plate_solve round-trip, which needs an API key not set in this environment); real ZTF demo re-run post-fix, confirming identical baseline/ZOGY tracklet counts (133/667, both recovering the same 2 known objects) since the fix never touched detection positions. commit - 9bb25233085f9ed8c126121aec1558e40f34edbe commit + dd2afcad8fb31f92714d78933a8e2aa88a22de0a blob - 720949b328aa480c1ef28bae6e4e4756171c740a blob + 21a81c867be7bc24fc6b742b60d8623c241535e1 --- INVESTIGATION_LOG.md +++ INVESTIGATION_LOG.md @@ -154,3 +154,86 @@ existing WCS ignored) solved successfully: recovered service; a retry succeeded, so `plate_solve` callers should be prepared to retry on transient server errors rather than assume a single failed upload means the service is unusable. + +## `detect_sources`'s flux has been silently wrong since the beginning: a transposed aperture + +Building `find_variable_sources` required trusting `detect_sources`'s +`flux` for the first time in a *quantitative* way — every earlier use +(`link_candidates`, tracklet building, WCS calibration, cross-matching) +only depends on its `x`/`y` positions. That new dependency surfaced a bug +that had been present, silently, since `detect_sources` was first written: +its measured flux was wrong by orders of magnitude for almost any real +source. + +Root cause: `Photometry.jl` (v0.9.8) is internally inconsistent between +its own two pieces. `PeakMesh`'s `extract_sources` reports positions in +the standard Cartesian sense — `to_nt(ci) = (x=ci[2], y=ci[1], ...)`, i.e. +`x` is the array's *second* dimension (column), `y` its *first* (row). +`CircularAperture`, in the same package, does the opposite internally — +`bounds`/`overlap` treat its own `.x` field as indexing the array's +*first* dimension and `.y` the second. `detect_sources` built +`CircularAperture(row.x, row.y, aperture_radius)` directly from +`PeakMesh`'s output, so every aperture was centred at the transposed +pixel — correct only when a source's row and column indices happened to +coincide (the diagonal), or invisibly wrong on a square canvas at a +generic position, or (on a non-square canvas, or once truly out of the +transposed array's bounds) landing on pure background instead. + +Found by comparing `detect_sources`'s measured flux against the closed-form +enclosed-energy integral for an isolated, noise-free synthetic Gaussian on +a non-square canvas at an asymmetric position: measured flux came back +~140x too small. Swapping the aperture's constructor arguments +(`CircularAperture(row.y, row.x, aperture_radius)`) brought it to within +~1.4% of the analytic value — geometric quantization error, not a +remaining bug. `light_curve` (`src/rotation.jl`) was checked the same way +and is *not* affected: it deliberately never permutes its image (see its +own docstring), and `WCS.world_to_pix`'s returned `(x, y)` already matches +raw FITS `(NAXIS1, NAXIS2)` order — which happens to be exactly what +`CircularAperture` expects internally, by coincidence of two conventions +cancelling out rather than by any intentional match. + +This had zero effect on every result validated so far in this project — +`run_pipeline`'s real-data tracklet counts (133 baseline / 667 ZOGY, both +recovering the same 2 known objects) depend only on detected *positions*, +never on `flux` — but it fully invalidated the first real-data +measurements made *while building* `find_variable_sources` earlier the +same session, before this was found: a claimed 29.4% median +`flux_err/flux` (actual, post-fix: **2.56%**), a claimed "only 2 of 120 +stars clear a 10% S/N floor" (actual: **119 of 119**), and a claimed +8-40% ensemble-vs-`MAGZP` mismatch explained by low-S/N selection bias +(the mismatch is real, but stable at 9-13% with or without an S/N cut — +not a selection effect at all; see below). Every constant and docstring +claim built on those numbers was rewritten against the corrected +measurements rather than left standing. + +## What was actually driving the ~10% photometric-scale mismatch, once flux was measured correctly + +With flux measured correctly, `photometric_scale`'s ensemble ratio against +each frame's own `MAGZP` zeropoint still disagreed by 9-13% — but now +*independent* of the S/N cut, ruling out selection bias as the cause. +Checked against each frame's `SEEING` header value (1.805-2.009 px across +the 5 real frames): the frames with better (smaller) seeing than frame 1 +all showed `photometric_scale < 1`, consistent with the standard "aperture +correction" effect — a fixed-radius aperture (`aperture_radius`, default 3 +px, comparable to ZTF's own seeing) encloses a larger fraction of a star's +total flux when the PSF is more concentrated. `photometric_scale`'s +ensemble-differential approach corrects for this automatically (it only +needs frame-to-frame consistency, not a causal model), so no code change +was needed here — only the docstring's explanation, which had cited the +now-debunked selection-bias story. + +## `find_variable_sources`'s chi-squared test has a real, measured false-positive floor on real data + +Even with correct flux and photometric normalization, real ZTF data's +reduced-chi2 distribution over 119 matched stationary stars has a heavy +tail: 23% exceed a threshold of 3 (the textbook-reasonable default), +13% still exceed 20. The likely cause: `detect_sources` positions each +detection at `PeakMesh`'s integer peak pixel, not a sub-pixel centroid, so +a 1-pixel jitter between frames — from noise, or a slightly different PSF +realization — against a small aperture (3 px, close to the PSF's own +scale) produces a real, but spurious, flux swing frame to frame. Raising +`chi2_threshold`'s default to 50 cuts the false-positive rate to 8% on +this dataset — better, but not clean, and documented as such in +`find_variable_sources`'s own docstring rather than presented as solved. +The actual fix (forced, sub-pixel-centroided photometry instead of +peak-position aperture photometry) is a larger change, not attempted here. blob - 37b103dafba9f9f78a2a42438c7918ed11920ae3 blob + c2cd7f2d371d9490ac8ba54a4d61056171942ed6 --- README.md +++ README.md @@ -69,9 +69,16 @@ helps. `find_variable_sources`/`search_field` (stationary, flux-varying source detection) and `fit_moffat_psf` (`estimate_psf`'s analytic-PSF fallback) -are implemented and validated only against synthetic data so far — not -yet exercised against real survey frames the way the rest of this list -has been. +are implemented, tested against synthetic data, and — for +`find_variable_sources`'s photometric normalization, S/N floor, and +`chi2_threshold` default, and for `fit_moffat_psf`'s recovered PSF +width — calibrated directly against real ZTF data (see +[`INVESTIGATION_LOG.md`](INVESTIGATION_LOG.md), including a real bug in +`detect_sources`'s own flux measurement this calibration work found and +fixed). Not yet exercised end to end, via `search_field`, against a real +field with an independently-confirmed variable star as ground truth — the +one real dataset checked so far has only one catalogued (VSX) variable in +its footprint, too faint to serve as a useful positive control. On that real dataset (field 451, 2019-10-23), the undifferenced baseline finds 133 tracklets and recovers both known objects in the field (2002 @@ -109,6 +116,21 @@ so far, all fixed with regression tests) and how each `length(fits_paths) - 1`) when using the ZOGY path, or no tracklet will ever be reachable if any frame gets gated — `examples/real_data_demo.jl` does this. +- **`find_variable_sources` has a real, measured false-positive floor on + real single-epoch aperture photometry.** Peak-pixel (not sub-pixel + centroid) positions mean a 1-pixel jitter against a small aperture can + look like genuine variability; on real ZTF data even a generous + `chi2_threshold` still flags several times more stars than the true + stellar variable fraction (see `find_variable_sources`'s docstring and + [`INVESTIGATION_LOG.md`](INVESTIGATION_LOG.md) for the measured rate). + Treat a candidate as needing independent confirmation (a catalog match + or a recovered period), not as self-evidently real. +- **`crossmatch_catalog(...; :vsx)`/`(...; :simbad)` query one candidate + at a time.** Migrated off the CDS X-Match service (extended, total + outages — see [`INVESTIGATION_LOG.md`](INVESTIGATION_LOG.md)) to direct + SIMBAD/VizieR TAP queries, which don't offer X-Match's single-batched-request + shape; a large candidate list means that many requests. `:skybot` is + unaffected (a different service, always queried this way). ## Example: real data blob - 75c81c5b485a170e68fd8f08dffddca42b8bfb0f blob + 2c9143232872c2a7de4adfba3d22e9f7ce278a71 --- src/AsteroidPipeline.jl +++ src/AsteroidPipeline.jl @@ -30,6 +30,6 @@ include("rotation.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 + search_field, find_variable_sources, variability_chi2, photometric_scale end # module AsteroidPipeline blob - e190aebcf5056e78c0ede4589b2d9654297aef77 blob + 0f489f8392a022fd9bc05c6626a77f7ebc863844 --- src/crossmatch.jl +++ src/crossmatch.jl @@ -1,5 +1,5 @@ -const _CDS_XMATCH_URL = "https://cdsxmatch.u-strasbg.fr/xmatch/api/v1/sync" -const _CDS_CATALOG_NAMES = Dict(:vsx => "vizier:B/vsx/vsx", :simbad => "simbad") +const _SIMBAD_TAP_URL = "https://simbad.cds.unistra.fr/simbad/sim-tap/sync" +const _VIZIER_TAP_URL = "https://tapvizier.cds.unistra.fr/TAPVizieR/tap/sync" const _SKYBOT_URL = "https://vo.imcce.fr/webservices/skybot/skybotconesearch_query.php" @@ -14,48 +14,117 @@ known objects from candidates warranting human verific J2000), such as the table returned by [`astrometric_calibrate`](@ref). SkyBoT additionally requires an `epoch` field (Julian Date) on each row, since it computes solar-system-object ephemerides for a specific instant -rather than querying a static catalog; `:vsx` and `:simbad` are queried in -a single batched request via the CDS X-Match service. +rather than querying a static catalog. `:vsx` and `:simbad` are each +queried directly against their own CDS TAP service (ADQL cone search), +**one request per candidate** — not the single batched request the CDS +X-Match service used to offer, since X-Match itself became unreliable +(extended, total outages; see `INVESTIGATION_LOG.md`) while the +underlying SIMBAD and VizieR TAP services stayed up. The trade-off is real +(N candidates means N requests, not 1) and matters for a large candidate +list — batch by querying a shared sky region directly via TAP if that +becomes a bottleneck; not done here since it wasn't yet. Returns a table with one row per (candidate, catalog match) pair found -within `radius`. Candidate `id`s absent from the returned table have no -known counterpart in `catalog`. +within `radius`; `:vsx` and `:simbad` return different columns (VSX +carries variability class/magnitude/period, SIMBAD doesn't), matching +what each catalog actually offers, but both always include `id`, `ra`, +`dec`, and `distance_arcsec`. Candidate `id`s absent from the returned +table have no known counterpart in `catalog`. """ function crossmatch_catalog(candidates, catalog::Symbol; radius::Real) if catalog === :skybot return _crossmatch_skybot(candidates, radius) - elseif haskey(_CDS_CATALOG_NAMES, catalog) - return _crossmatch_cds(candidates, _CDS_CATALOG_NAMES[catalog], radius) + elseif catalog === :simbad + return _crossmatch_simbad(candidates, radius) + elseif catalog === :vsx + return _crossmatch_vsx(candidates, radius) else throw(ArgumentError("unknown catalog $(repr(catalog)); expected :skybot, :vsx, or :simbad")) end end -function _crossmatch_cds(candidates, cds_name::AbstractString, radius::Real) - 0 < radius <= 180 || - throw(ArgumentError("radius must be in (0, 180] arcsec for the CDS X-Match service")) +""" + _tap_query(url, query) -> CSV.File - upload = IOBuffer() - println(upload, "id,ra,dec") - for c in candidates - println(upload, "$(c.id),$(c.ra),$(c.dec)") - end - seekstart(upload) - +Run `query` (an ADQL string) as a synchronous TAP query against `url`, +returning the CSV response parsed by `CSV.File`. Shared by +[`_crossmatch_simbad`](@ref) and [`_crossmatch_vsx`](@ref). +""" +function _tap_query(url::AbstractString, query::AbstractString) body = HTTP.Form(Dict( - "request" => "xmatch", - "distMaxArcsec" => string(radius), - "RESPONSEFORMAT" => "csv", - "cat1" => HTTP.Multipart("candidates.csv", upload, "text/csv"), - "colRA1" => "ra", - "colDec1" => "dec", - "cat2" => cds_name, - )) - response = HTTP.post(_CDS_XMATCH_URL, [], body) + "REQUEST" => "doQuery", "LANG" => "ADQL", "FORMAT" => "csv", "QUERY" => query)) + response = HTTP.post(url, [], body) + return CSV.File(response.body) +end - return Table(CSV.File(response.body)) +""" + _cds_cone_query(select, table, ra_col, dec_col, ra, dec, radius_deg) -> String + +An ADQL synchronous cone-search query: rows of `table` within +`radius_deg` of `(ra, dec)`, plus a `distance_arcsec` column (via ADQL's +`DISTANCE`, in degrees, converted here) — shared by +[`_crossmatch_simbad`](@ref) and [`_crossmatch_vsx`](@ref), which differ +only in `select`/`table`/coordinate column names. +""" +function _cds_cone_query(select::AbstractString, table::AbstractString, + ra_col::AbstractString, dec_col::AbstractString, + ra::Real, dec::Real, radius_deg::Real) + point = "POINT('ICRS', $ra_col, $dec_col)" + return "SELECT $select, DISTANCE($point, POINT('ICRS', $ra, $dec)) * 3600 AS distance_arcsec " * + "FROM $table WHERE CONTAINS($point, CIRCLE('ICRS', $ra, $dec, $radius_deg)) = 1" end +_nan_if_missing(x) = x === missing ? NaN : Float64(x) + +function _crossmatch_simbad(candidates, radius::Real) + 0 < radius <= 180 || throw(ArgumentError("radius must be in (0, 180] arcsec")) + radius_deg = radius / 3600 + + id = Int[]; name = String[]; ra = Float64[]; dec = Float64[]; distance_arcsec = Float64[] + for c in candidates + query = _cds_cone_query("main_id, ra, dec", "basic", "ra", "dec", c.ra, c.dec, radius_deg) + for row in _tap_query(_SIMBAD_TAP_URL, query) + push!(id, c.id) + push!(name, String(row.main_id)) + push!(ra, row.ra) + push!(dec, row.dec) + push!(distance_arcsec, row.distance_arcsec) + end + end + + return Table(; id, name, ra, dec, distance_arcsec) +end + +function _crossmatch_vsx(candidates, radius::Real) + 0 < radius <= 180 || throw(ArgumentError("radius must be in (0, 180] arcsec")) + radius_deg = radius / 3600 + + id = Int[]; name = String[]; class = String[]; ra = Float64[]; dec = Float64[] + mag_max = Float64[]; mag_min = Float64[]; period = Float64[]; distance_arcsec = Float64[] + for c in candidates + query = _cds_cone_query("Name, Type, RAJ2000, DEJ2000, max, min, Period", "\"B/vsx/vsx\"", + "RAJ2000", "DEJ2000", c.ra, c.dec, radius_deg) + for row in _tap_query(_VIZIER_TAP_URL, query) + push!(id, c.id) + # VizieR's TAP service returns Name/Type as fixed-width, + # space-padded strings (e.g. "RS" arrives as "RS" followed by + # 28 spaces) — found via a real query, not documented anywhere + # obvious; strip or every string comparison against these + # silently fails. + push!(name, String(strip(row.Name))) + push!(class, String(strip(row.Type))) + push!(ra, row.RAJ2000) + push!(dec, row.DEJ2000) + push!(mag_max, _nan_if_missing(row.max)) + push!(mag_min, _nan_if_missing(row.min)) + push!(period, _nan_if_missing(row.Period)) + push!(distance_arcsec, row.distance_arcsec) + end + end + + return Table(; id, name, class, ra, dec, mag_max, mag_min, period, distance_arcsec) +end + function _crossmatch_skybot(candidates, radius::Real) radius_deg = radius / 3600 blob - 2e302e53c178203fbb3d4fc01e332dcc3106d494 blob + 5481c545ad29ec7a14f867eaf97b348f077de8e7 --- src/detection.jl +++ src/detection.jl @@ -1,6 +1,7 @@ """ detect_sources(image::AbstractMatrix{<:Real}; threshold::Real, - box_size::NTuple{2,Integer}=(5, 5), aperture_radius::Real=3.0) + box_size::NTuple{2,Integer}=(5, 5), aperture_radius::Real=3.0, + gain::Union{Nothing,Real}=nothing) Detect point sources in a FITS frame. @@ -13,15 +14,26 @@ detection's flux is then measured with circular apertu radius `aperture_radius` pixels on the background-subtracted image. Returns a table with columns `x`, `y` (pixel position), `peak` (background- -subtracted peak pixel value), `flux` (aperture sum), and `flux_err` -(`Photometry.photometry`'s propagated aperture error from the uniform -per-pixel `noise` used for detection — not a full per-pixel variance map, -the same caveat [`light_curve`](@ref) documents), to be consumed by -[`link_candidates`](@ref) for inter-frame motion matching or -[`find_variable_sources`](@ref) for flux-variability matching. +subtracted peak pixel value), `flux` (aperture sum), and `flux_err`. By +default (`gain=nothing`) `flux_err` is `Photometry.photometry`'s propagated +aperture error from the uniform per-pixel `noise` used for detection alone +— not a full per-pixel variance map, the same caveat [`light_curve`](@ref) +documents. Passing `gain` (electrons/ADU, from the frame's own `GAIN` +header keyword) adds each pixel's own Poisson (shot) noise, +`sqrt(noise^2 + max(pixel, 0) / gain)`: background noise alone +underestimates a bright star's real flux uncertainty — on real ZTF data +this barely moves the *median* `flux_err/flux` (2.97% vs 3.16%, most +detections being near-threshold and background-noise-dominated either +way), but for the single brightest star in that same frame, shot noise +was ~10x the background-only estimate (0.039% vs 0.004%) — exactly where +underestimating the error bar would most distort +[`find_variable_sources`](@ref)'s chi-squared test. Left `nothing` +(background noise only) for [`link_candidates`](@ref), which only uses +`x`/`y` and has no use for a flux uncertainty at all. """ function detect_sources(image::AbstractMatrix{<:Real}; threshold::Real, - box_size::NTuple{2,<:Integer}=(5, 5), aperture_radius::Real=3.0) + box_size::NTuple{2,<:Integer}=(5, 5), aperture_radius::Real=3.0, + gain::Union{Nothing,Real}=nothing) background, noise = estimate_background(image; location=SourceExtractorBackground(), rms=MADStdRMS()) subtracted = image .- background @@ -30,11 +42,26 @@ function detect_sources(image::AbstractMatrix{<:Real}; end finder = PeakMesh(box_size, threshold) - error_map = fill(noise, size(image)) - peaks = extract_sources(finder, subtracted, error_map) + detection_error_map = fill(noise, size(image)) + peaks = extract_sources(finder, subtracted, detection_error_map) - apertures = [CircularAperture(row.x, row.y, aperture_radius) for row in peaks] - photom = isempty(apertures) ? nothing : photometry(apertures, subtracted, error_map) + # `PeakMesh` reports x/y in the standard Cartesian sense (x=column, + # i.e. the array's 2nd dimension; y=row, the 1st — see + # Photometry.jl's own `extract_sources`, `to_nt(ci) = (x=ci[2], + # y=ci[1], ...)`). `CircularAperture`, in the very same package, + # does the opposite internally (its `x` field indexes the array's + # *1st* dimension, `y` the 2nd — see `bounds`/`overlap` in + # Photometry.jl's circular.jl). Passing `(row.x, row.y)` straight + # through silently centers the aperture at the transposed pixel + # whenever the true position isn't on the row==column diagonal — + # confirmed by comparing a measured flux against the analytic + # enclosed-energy integral for an isolated Gaussian, which came back + # ~140x too small before this swap and matched to ~1% after it. + # `(row.y, row.x)` here is the fix, not a second bug. + apertures = [CircularAperture(row.y, row.x, aperture_radius) for row in peaks] + photom_error_map = gain === nothing ? detection_error_map : + sqrt.(noise^2 .+ max.(subtracted, 0.0) ./ gain) + photom = isempty(apertures) ? nothing : photometry(apertures, subtracted, photom_error_map) flux = isempty(apertures) ? Float64[] : [row.aperture_sum for row in photom] flux_err = isempty(apertures) ? Float64[] : [row.aperture_sum_err for row in photom] blob - 44855504a9f86db7087ccfe19884bf2b5a801e40 blob + 320b8e8626a390617bfaec3e4baafa730d793a7a --- src/pipeline.jl +++ src/pipeline.jl @@ -5,7 +5,8 @@ aperture_radius::Real=3.0, max_speed::Real=Inf, match_radius::Real=2.0, min_frames::Integer=length(fits_paths), reference=nothing, psf_threshold::Real=20.0, psf_min_separation::Real=40.0, - quality_max_std::Real=1.5, plate_solve_api_key::Union{Nothing,AbstractString}=nothing) + quality_max_std::Real=1.5, plate_solve_api_key::Union{Nothing,AbstractString}=nothing, + photometric_outlier_threshold::Real=0.2) Run the local stages of the pipeline — detection, linking, and astrometric calibration — on a time-ordered sequence of FITS frames from @@ -65,6 +66,24 @@ actually separates them here. The default of `1.5` sit both sides of that real gap; it is calibrated from that one dataset, not a universal constant — tune per survey/conditions. +The raw (no-`reference`) path has no equivalent gate built into +`S_corr`'s own statistics, but a related signal is available there: +[`photometric_scale`](@ref)'s per-frame flux-scale factor (computed after +detection, from every frame's own detections). If any frame's factor +relative to the field's median deviates by more than +`photometric_outlier_threshold` (default `0.2`, i.e. 20%), a warning is +emitted — this is only a warning, never an automatic exclusion, unlike +`quality_max_std`, precisely to avoid repeating that gate's own +combinatorial side effect on `link_candidates`'s `min_frames` (see +`INVESTIGATION_LOG.md`). Unlike `quality_max_std`, this has **not** been +validated against a real, independently-confirmed anomaly: on the same 5 +real ZTF frames used to calibrate `quality_max_std` (one of which is a +confirmed likely passing cloud), every frame's raw-path photometric scale +stayed within ~3-13% of the field median — well under the default +threshold, so this signal did not (and, on this dataset, could not have) +flagged that frame. It is included as a plausible general-purpose check, +not a proven one; treat the default threshold as unvalidated until it is. + If a frame's header has no WCS, [`load_wcs`](@ref) raises an error; passing `plate_solve_api_key` (a nova.astrometry.net API key) makes that frame fall back to [`plate_solve`](@ref) instead of failing outright — @@ -88,10 +107,12 @@ function run_pipeline(fits_paths::AbstractVector{<:Abs aperture_radius::Real=3.0, max_speed::Real=Inf, match_radius::Real=2.0, min_frames::Integer=length(fits_paths), reference=nothing, psf_threshold::Real=20.0, psf_min_separation::Real=40.0, - quality_max_std::Real=1.5, plate_solve_api_key::Union{Nothing,AbstractString}=nothing) + quality_max_std::Real=1.5, plate_solve_api_key::Union{Nothing,AbstractString}=nothing, + photometric_outlier_threshold::Real=0.2) detections_per_frame, wcs_per_frame, timestamps = _detect_all_frames( fits_paths; timestamp_key, threshold, box_size, aperture_radius, - reference, psf_threshold, psf_min_separation, quality_max_std, plate_solve_api_key) + reference, psf_threshold, psf_min_separation, quality_max_std, plate_solve_api_key, + photometric_outlier_threshold) tracklets = link_candidates(detections_per_frame, timestamps; max_speed, match_radius, min_frames) return astrometric_calibrate(tracklets, wcs_per_frame, timestamps) @@ -108,15 +129,23 @@ links these into movers) and [`search_field`](@ref) (w looks for stationary variable sources) — factored out so both share one detection pass over `fits_paths` rather than repeating the expensive reprojection/PSF/ZOGY work. See `run_pipeline`'s docstring for the -meaning of every keyword; behaviour here is identical to what -`run_pipeline` did inline before this split. +meaning of every keyword. On the raw (no-`reference`) path, each frame's +own `GAIN` header keyword (default `1.0` if absent) is passed to +[`detect_sources`](@ref) so `flux_err` includes source shot noise, not +just background noise — needed for [`find_variable_sources`](@ref)'s +chi-squared test to have a realistic error bar to test against. Not +applied on the `reference` (ZOGY) path, since `S_corr` there is already a +normalized detection-significance map, not physical counts. On that same +raw path, `photometric_outlier_threshold` is forwarded to a post-loop +[`photometric_scale`](@ref) check — see `run_pipeline`'s docstring. """ function _detect_all_frames(fits_paths::AbstractVector{<:AbstractString}; timestamp_key::AbstractString="MJD-OBS", threshold::Real=5.0, box_size::NTuple{2,<:Integer}=(5, 5), aperture_radius::Real=3.0, reference=nothing, psf_threshold::Real=20.0, psf_min_separation::Real=40.0, - quality_max_std::Real=1.5, plate_solve_api_key::Union{Nothing,AbstractString}=nothing) + quality_max_std::Real=1.5, plate_solve_api_key::Union{Nothing,AbstractString}=nothing, + photometric_outlier_threshold::Real=0.2) detections_per_frame = [] wcs_per_frame = WCSTransform[] timestamps = Float64[] @@ -185,13 +214,29 @@ function _detect_all_frames(fits_paths::AbstractVector end push!(wcs_per_frame, reference.wcs) end - push!(detections_per_frame, detect_sources(image; threshold, box_size, aperture_radius)) + # gain only applies to the raw path: S_corr (the reference + # path) is already a normalized detection-significance map, + # not physical counts, so Poisson noise doesn't apply to it. + push!(detections_per_frame, + detect_sources(image; threshold, box_size, aperture_radius, + gain=(reference === nothing ? gain : nothing))) mjd = read_key(hdu, timestamp_key)[1] push!(timestamps, mjd + 2400000.5) end end + if reference === nothing && length(fits_paths) >= 2 + scales = photometric_scale(detections_per_frame) + reference_scale = median(scales) + for (k, path) in enumerate(fits_paths) + deviation = abs(scales[k] - reference_scale) / reference_scale + if deviation > photometric_outlier_threshold + @warn "frame's photometric scale deviates from the field's median" path scale=scales[k] reference_scale deviation photometric_outlier_threshold + end + end + end + return detections_per_frame, wcs_per_frame, timestamps end @@ -199,7 +244,7 @@ end search_field(fits_paths; , variability_position_tolerance::Real=2.0, variability_min_frames::Integer=length(fits_paths), - variability_chi2_threshold::Real=3.0) + variability_chi2_threshold::Real=50.0) -> (movers=, variables=
) Run [`run_pipeline`](@ref)'s asteroid-candidate search and @@ -211,11 +256,16 @@ the efficient choice for a real observing run that wan since `run_pipeline` alone would need a second full pass over the same frames to also find variables. -All keywords through `plate_solve_api_key` are exactly `run_pipeline`'s -(see its docstring); `variability_position_tolerance`, -`variability_min_frames`, and `variability_chi2_threshold` are forwarded -to `find_variable_sources` as its `position_tolerance`, `min_frames`, and -`chi2_threshold`. +All keywords through `photometric_outlier_threshold` are exactly +`run_pipeline`'s (see its docstring); `variability_position_tolerance`, +`variability_min_frames`, `variability_chi2_threshold`, +`variability_normalize`, and `variability_max_relative_error` are +forwarded to [`find_variable_sources`](@ref) as its `position_tolerance`, +`min_frames`, `chi2_threshold`, `normalize`, and `max_relative_error`. +`variability_normalize` defaults to `reference === nothing` — off on the +ZOGY path, per `find_variable_sources`'s own docstring (`S_corr` isn't on +a physical flux scale, so ensemble photometric normalization doesn't +apply there). Returns a named tuple `(movers=..., variables=...)`, each an `astrometric_calibrate` candidate table (columns `id`, `frame`, `x`, `y`, @@ -231,12 +281,16 @@ function search_field(fits_paths::AbstractVector{<:Abs match_radius::Real=2.0, min_frames::Integer=length(fits_paths), reference=nothing, psf_threshold::Real=20.0, psf_min_separation::Real=40.0, quality_max_std::Real=1.5, plate_solve_api_key::Union{Nothing,AbstractString}=nothing, + photometric_outlier_threshold::Real=0.2, variability_position_tolerance::Real=2.0, variability_min_frames::Integer=length(fits_paths), - variability_chi2_threshold::Real=3.0) + variability_chi2_threshold::Real=50.0, + variability_normalize::Bool=(reference === nothing), + variability_max_relative_error::Real=0.10) detections_per_frame, wcs_per_frame, timestamps = _detect_all_frames( fits_paths; timestamp_key, threshold, box_size, aperture_radius, - reference, psf_threshold, psf_min_separation, quality_max_std, plate_solve_api_key) + reference, psf_threshold, psf_min_separation, quality_max_std, plate_solve_api_key, + photometric_outlier_threshold) tracklets = link_candidates(detections_per_frame, timestamps; max_speed, match_radius, min_frames) movers = astrometric_calibrate(tracklets, wcs_per_frame, timestamps) @@ -245,7 +299,9 @@ function search_field(fits_paths::AbstractVector{<:Abs detections_per_frame, timestamps; position_tolerance=variability_position_tolerance, min_frames=variability_min_frames, - chi2_threshold=variability_chi2_threshold) + chi2_threshold=variability_chi2_threshold, + normalize=variability_normalize, + max_relative_error=variability_max_relative_error) variables = astrometric_calibrate(variable_groups, wcs_per_frame, timestamps) return (movers=movers, variables=variables) blob - 331b15fbab956473fab4da0314ec1761a29c059a blob + ea0c2cf648312a3036dfb165a262541856ea8d67 --- src/variables.jl +++ src/variables.jl @@ -21,10 +21,86 @@ function variability_chi2(flux::AbstractVector{<:Real} end """ + _filter_high_snr(detections, max_relative_error) -> Vector + +`detections` (a table with `x`, `y`, `flux`, `flux_err` columns) restricted +to rows with positive flux and `flux_err / flux < max_relative_error` — the +S/N floor shared by [`photometric_scale`](@ref) and +[`find_variable_sources`](@ref), factored out so both apply exactly the +same cut. +""" +_filter_high_snr(detections, max_relative_error::Real) = + [row for row in detections if row.flux > 0 && row.flux_err / row.flux < max_relative_error] + +""" + photometric_scale(detections_per_frame; position_tolerance::Real=2.0, + max_relative_error::Real=0.10, min_stars::Integer=5) -> Vector{Float64} + +Per-frame photometric scale factor, relative to frame 1, by ensemble +differential photometry — the survey-agnostic way to put every frame's +flux on a common scale without depending on any zeropoint header keyword +(real IASC campaign headers are unverified — see `README.md`). + +For each frame `k`, matches stars against frame 1 by position (within +`position_tolerance` pixels, via [`link_candidates`](@ref)'s +`_closest_detection`), and takes the **median** of `flux_1 / flux_k` over +those matches — restricted, in both frames, to stars with +`flux_err / flux < max_relative_error`, so one noisy or blended star can't +swing the whole frame's scale. + +Measured directly on real ZTF data (field 451, 2019-10-23, 119 matched +stars): the ensemble ratio disagreed with the frames' own `MAGZP` +zeropoints by 9-13%, essentially *unchanged* by the S/N cut above +(-8.9% to -13.2% either way) — so this is not a low-S/N selection effect. +It tracks each frame's own `SEEING` instead (1.805-2.009 px across the +five frames): a fixed-radius aperture (`aperture_radius` in +[`detect_sources`](@ref)) encloses a seeing-dependent fraction of a star's +total PSF flux, the standard "aperture correction" effect in photometry. +`photometric_scale`'s ensemble approach corrects for exactly this kind of +uniform per-frame multiplicative offset, whatever its cause — it does not +need to know *why* frame `k`'s stars all read low or high, only that they +do, consistently. + +If fewer than `min_stars` matches survive the S/N cut for some frame, +that frame's factor is left at `1.0` (no correction) with a `@warn`, +rather than trusting a scale derived from a handful of stars. + +Returns a vector the same length as `detections_per_frame`; index 1 is +always `1.0` by construction (frame 1 is its own reference). +""" +function photometric_scale(detections_per_frame; position_tolerance::Real=2.0, + max_relative_error::Real=0.10, min_stars::Integer=5) + nframes = length(detections_per_frame) + scales = ones(Float64, nframes) + nframes < 2 && return scales + + ref = _filter_high_snr(detections_per_frame[1], max_relative_error) + for k in 2:nframes + candidates = _filter_high_snr(detections_per_frame[k], max_relative_error) + ratios = Float64[] + for d1 in ref + best = _closest_detection(candidates, d1.x, d1.y, position_tolerance) + best === nothing && continue + push!(ratios, d1.flux / best.flux) + end + + if length(ratios) < min_stars + @warn "too few high-S/N matched stars to measure frame's photometric scale; leaving uncorrected" frame=k n_matches=length(ratios) min_stars max_relative_error + continue + end + scales[k] = median(ratios) + end + + return scales +end + +""" find_variable_sources(detections_per_frame, timestamps; position_tolerance::Real=2.0, min_frames::Integer=length(detections_per_frame), - chi2_threshold::Real=3.0) + chi2_threshold::Real=50.0, + normalize::Bool=true, + max_relative_error::Real=0.10) Match source detections across frames by consistent *position* (zero assumed motion) rather than [`link_candidates`](@ref)'s linear-motion @@ -38,39 +114,75 @@ matching [`link_candidates`](@ref)'s convention (used frame order meaningful for periodicity follow-up, not for any position prediction). -For each detection in frame 1, the closest detection within -`position_tolerance` pixels of that *same* position in every other frame -is matched (via [`link_candidates`](@ref)'s `_closest_detection` — no -velocity model, since a variable star does not move between frames). -Groups reaching at least `min_frames` matched points are then tested -against [`variability_chi2`](@ref)'s constant-flux null hypothesis; a -group is kept only if its **reduced** chi-squared (`chi2 / dof`) exceeds -`chi2_threshold`. +Two corrections are applied before any variability test: -What "keep only if variable" filters depends on how `detections_per_frame` -was produced: +- **Low-S/N floor** (`max_relative_error`, default 10%): detections with + `flux_err / flux` at or above this are dropped before matching, in every + frame — a safety net against reporting "variability" from measurements + too imprecise to support the claim. +- **Photometric normalization** (`normalize`, default `true`): a uniform + per-frame flux-scale offset (see [`photometric_scale`](@ref) — on real + ZTF data this tracked each frame's own seeing, a fixed-aperture-radius + effect, not a code defect) is corrected before the chi-squared test, so + a constant star doesn't read as variable just because one exposure's + aperture happened to enclose a different fraction of its PSF. Set + `normalize=false` to skip this (e.g. if `flux` is already calibrated). +For each detection in frame 1 (after the S/N floor), the closest detection +within `position_tolerance` pixels of that *same* position in every other +frame (after its own S/N floor) is matched. Groups reaching at least +`min_frames` matched points are tested against [`variability_chi2`](@ref)'s +constant-flux null hypothesis; a group is kept only if its **reduced** +chi-squared (`chi2 / dof`) exceeds `chi2_threshold`. + +`chi2_threshold`'s default is deliberately not close to the naive "3" a +formally-correct chi-squared test would suggest. Measured on real ZTF +field 451 (119 matched, high-S/N stationary stars, after normalization): +the reduced-chi2 distribution has a real, heavy tail — 23% of ordinary +stars exceed a threshold of 3, and even a threshold of 20 still flags +13% — almost certainly because `detect_sources` finds a peak-detected +integer pixel each frame, not a sub-pixel centroid, so a 1 px jitter +against a small `aperture_radius` (default 3 px, comparable to a typical +seeing FWHM) shows up as a real, if spurious, flux change. `50.0` cuts +that to 8% on the same data — still elevated above the true stellar +variable fraction, not a clean cut, and honestly reported as such rather +than tuned to look better than it is: treat any candidate here as +requiring independent confirmation (a VSX/SIMBAD match via +[`crossmatch_catalog`](@ref), or a period recovered by +[`recover_rotation_period`](@ref)), not as self-evidently real. Properly +fixing the underlying cause — forced, sub-pixel-centroided photometry +instead of peak-position aperture photometry — is future work, not +attempted here. + +With the default `min_frames` (every frame must match at high S/N), a +field where few sources clear the S/N floor in every frame will correctly +return few or no candidates rather than a list built on unreliable +photometry — lower `max_relative_error` or `min_frames` deliberately if a +shallower search is wanted, rather than reading an empty result as "no +code path found here". + +What the chi2 test is discriminating against depends on how +`detections_per_frame` was produced: + - **ZOGY path** (`detections_per_frame` from `S_corr` — `run_pipeline` or `search_field` given a `reference`): a constant star differences to zero and is never detected in `S_corr` at all, so a matched stationary group is already a variable/transient by construction; the chi2 test here is a second filter against the known subtraction-residual failure mode at bright stars (see `INVESTIGATION_LOG.md`), not the primary - discriminant. + discriminant. `normalize`/`photometric_scale` do not apply meaningfully + here, since `S_corr` isn't on a physical flux scale to begin with — + pass `normalize=false` on this path. - **Raw path** (no `reference`): every star above `threshold` is detected - in every frame, so the chi2 test is the only thing separating a genuine - variable from an ordinary constant star. + in every frame, so the chi2 test (with normalization and the S/N floor) + is the only thing separating a genuine variable from an ordinary + constant star. -`chi2_threshold`'s right value is data-dependent (like `threshold` and -`quality_max_std` elsewhere in this package) — it assumes `flux_err` is -the propagated-aperture-error estimate `detect_sources` provides (from a -uniform per-pixel noise, not a full variance map), so treat the absolute -chi2 scale as approximate and tune per survey. - Returns a vector of matched-point groups, each a vector of `(frame, x, y, flux, flux_err)` named tuples sorted by frame — the same -shape [`link_candidates`](@ref) returns (with two extra fields), so a -group can be passed directly to [`astrometric_calibrate`](@ref). +shape [`link_candidates`](@ref) returns (with two extra fields, and with +`flux`/`flux_err` already normalized when `normalize=true`), so a group +can be passed directly to [`astrometric_calibrate`](@ref). For periodicity, build `(times, flux)` from a returned group's own `frame`/`flux` fields and `timestamps`, and pass to @@ -79,7 +191,9 @@ For periodicity, build `(times, flux)` from a returned function find_variable_sources(detections_per_frame, timestamps; position_tolerance::Real=2.0, min_frames::Integer=length(detections_per_frame), - chi2_threshold::Real=3.0) + chi2_threshold::Real=50.0, + normalize::Bool=true, + max_relative_error::Real=0.10) nframes = length(detections_per_frame) nframes == length(timestamps) || throw(ArgumentError("detections_per_frame and timestamps must have the same length")) @@ -88,15 +202,20 @@ function find_variable_sources(detections_per_frame, t groups = Vector{Point}[] nframes < 1 && return groups - for d1 in detections_per_frame[1] + scales = normalize ? + photometric_scale(detections_per_frame; position_tolerance, max_relative_error) : + ones(Float64, nframes) + filtered = [_filter_high_snr(d, max_relative_error) for d in detections_per_frame] + + for d1 in filtered[1] group = Point[(frame=1, x=Float64(d1.x), y=Float64(d1.y), - flux=Float64(d1.flux), flux_err=Float64(d1.flux_err))] + flux=Float64(d1.flux) * scales[1], flux_err=Float64(d1.flux_err) * scales[1])] for k in 2:nframes - best = _closest_detection(detections_per_frame[k], d1.x, d1.y, position_tolerance) + best = _closest_detection(filtered[k], d1.x, d1.y, position_tolerance) best === nothing && continue push!(group, (frame=k, x=Float64(best.x), y=Float64(best.y), - flux=Float64(best.flux), flux_err=Float64(best.flux_err))) + flux=Float64(best.flux) * scales[k], flux_err=Float64(best.flux_err) * scales[k])) end length(group) < min_frames && continue blob - bdaaa5696ebca2069a9f947330a155fb60e619ff blob + 59a8d96fd87c5a38fe5881f1ba33b9c26c3a87f5 --- test/runtests.jl +++ test/runtests.jl @@ -28,6 +28,15 @@ using Reproject @test best.flux > 0 @test best.flux_err > 0 + # gain=nothing (default) uses background noise only; a finite gain + # adds the bright source's own Poisson (shot) noise, which must + # only ever widen flux_err, never narrow it, and only for a real + # source (a pure-noise frame's flux_err is unaffected since there's + # no positive signal to contribute shot noise). + with_gain = detect_sources(image; threshold=5.0, gain=2.0) + best_gain = with_gain[argmax(with_gain.peak)] + @test best_gain.flux_err > best.flux_err + noise_only = 100.0 .+ 5.0 .* randn(32, 32) @test length(detect_sources(noise_only; threshold=100.0)) == 0 end @@ -121,6 +130,52 @@ using Reproject @test all(==(1), calibrated.id) end + @testset "photometric_scale / normalized find_variable_sources" begin + # 5 constant stars, well above the default S/N floor, whose flux + # in frame 2 is uniformly 0.8x frame 1's (standing in for a real + # zeropoint/transparency mismatch between exposures — measured as + # large as this on real ZTF frames, see INVESTIGATION_LOG.md) plus + # one genuinely variable source whose flux does *not* follow that + # scaling; frame 3 returns to frame 1's scale. + star_x = [10.0, 20.0, 30.0, 40.0, 50.0] + star_y = copy(star_x) + frame1 = Table(x=vcat(star_x, 70.0), y=vcat(star_y, 70.0), + flux=vcat(fill(1000.0, 5), 2000.0), flux_err=vcat(fill(5.0, 5), 10.0)) + frame2 = Table(x=vcat(star_x, 70.0), y=vcat(star_y, 70.0), + flux=vcat(fill(800.0, 5), 2000.0), flux_err=vcat(fill(4.0, 5), 10.0)) + frame3 = Table(x=vcat(star_x, 70.0), y=vcat(star_y, 70.0), + flux=vcat(fill(1000.0, 5), 2000.0), flux_err=vcat(fill(5.0, 5), 10.0)) + timestamps = [0.0, 1.0, 2.0] + + scales = photometric_scale([frame1, frame2, frame3]) + @test scales[1] == 1.0 + @test scales[2]≈1.25 rtol=1e-6 # 1000/800 + @test scales[3]≈1.0 rtol=1e-6 + + groups = find_variable_sources([frame1, frame2, frame3], timestamps) + @test length(groups) == 1 + @test all(p.x == 70.0 && p.y == 70.0 for p in groups[1]) + + # without normalization, frame 2's uniform 0.8x mismatch alone can + # push a merely-rescaled constant star's chi2 over threshold too — + # normalization is what keeps that from happening above, so + # disabling it must not find *fewer* candidates here. + raw_groups = find_variable_sources([frame1, frame2, frame3], timestamps; normalize=false) + @test length(raw_groups) >= length(groups) + end + + @testset "photometric_scale min_stars guard" begin + # only 2 matched stars per frame — below the default min_stars=5 — + # so the scale must stay uncorrected (1.0) rather than trust a + # factor derived from 2 stars (measured on real data: only 2 of + # 120 detections cleared the default S/N floor in every frame). + frame1 = Table(x=[10.0, 20.0], y=[10.0, 20.0], flux=[1000.0, 1000.0], flux_err=[5.0, 5.0]) + frame2 = Table(x=[10.0, 20.0], y=[10.0, 20.0], flux=[800.0, 800.0], flux_err=[4.0, 4.0]) + + scales = @test_logs (:warn, r"too few") match_mode = :any photometric_scale([frame1, frame2]) + @test scales == [1.0, 1.0] + end + @testset "astrometry" begin wcs = WCSTransform(2; crpix=[500.0, 500.0], crval=[150.0, 20.0], cdelt=[-1 / 3600, 1 / 3600], ctype=["RA---TAN", "DEC--TAN"]) @@ -563,16 +618,29 @@ using Reproject @testset "search_field" begin mktempdir() do dir - nx, ny = 80, 60 + # Larger canvas than the "run_pipeline" testset above needs + # (80x60): a bright-enough variable star to clear the S/N floor + # below (amplitude up to 8000) contaminates + # BackgroundMeshes' background/RMS estimate at that canvas + # size — found empirically as ~175-450 spurious detections + # per frame at 80x60, dropping to a clean 2 once the source is + # a small enough fraction of the frame — the same + # source-fraction-of-image effect already documented in + # INVESTIGATION_LOG.md's quality-gate test (600x600 for the + # same reason). Real ZTF frames (~10^6 px) aren't at risk of + # this; a small synthetic frame with a very bright source is. + nx, ny = 300, 225 wcs = WCSTransform(2; crpix=[nx / 2, ny / 2], crval=[150.0, 20.0], cdelt=[-1 / 3600, 1 / 3600], ctype=["RA---TAN", "DEC--TAN"]) # true asteroid track (as in the "run_pipeline" testset above) x0, y0, dx, dy = 20.0, 45.0, 5.0, -3.0 # stationary variable star, well separated from the track, - # with a flux swing far beyond photon noise - xv, yv = 60.0, 15.0 - var_amps = [200.0, 900.0, 200.0] + # with a flux swing far beyond photon noise and flux_err/flux + # comfortably under find_variable_sources's default 10% S/N + # floor in every frame. + xv, yv = 220.0, 150.0 + var_amps = [3000.0, 8000.0, 3000.0] mjd0 = 60000.0 Random.seed!(3) @@ -692,6 +760,34 @@ using Reproject end try + # a real VSX variable (Gaia DR3 2501737571391422464, type RS) + # inside the real_data_demo.jl field, verified independently + # via a direct VizieR TAP query during the migration off the + # CDS X-Match service (see INVESTIGATION_LOG.md) — a positive + # control, not just an absence-of-error check. + star = [(id=1, ra=36.2344, dec=2.06997)] + matches = crossmatch_catalog(star, :vsx; radius=5.0) + @test any(==("RS"), matches.class) + @test matches[1].id == 1 + catch e + e isa HTTP.Exceptions.HTTPError || e isa Base.IOError || rethrow() + @test_skip "network unavailable" + end + + try + # near the celestial pole, tiny radius: regression test for + # empty CDS TAP results (mirrors the SkyBoT empty-result case + # below) — table construction must not crash on zero rows. + empty_candidates = [(id=1, ra=0.0, dec=89.9)] + matches = crossmatch_catalog(empty_candidates, :simbad; radius=1.0) + @test length(matches) == 0 + @test isempty(matches.id) + catch e + e isa HTTP.Exceptions.HTTPError || e isa Base.IOError || rethrow() + @test_skip "network unavailable" + end + + try # A known object (127319 "2002 JB99") ~134" from this position # at this epoch (verified independently against SkyBoT). This # is a positive control: it caught a real bug where the Julian