commit - 9bb25233085f9ed8c126121aec1558e40f34edbe
commit + dd2afcad8fb31f92714d78933a8e2aa88a22de0a
blob - 720949b328aa480c1ef28bae6e4e4756171c740a
blob + 21a81c867be7bc24fc6b742b60d8623c241535e1
--- INVESTIGATION_LOG.md
+++ INVESTIGATION_LOG.md
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
`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
`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
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
-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"
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
"""
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.
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
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
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
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 —
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)
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[]
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
search_field(fits_paths; <all run_pipeline keywords>,
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=<table>, variables=<table>)
Run [`run_pipeline`](@ref)'s asteroid-candidate search and
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`,
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)
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
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
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
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"))
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
@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
@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"])
@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)
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
