commit 8ff5933c85c68fcb18d5be2ee05c05d4d2c6df00 from: ale date: Mon Aug 17 04:13:52 2026 UTC Fix sub-pixel centroiding, auto min_frames, batched cross-matches; validate against real variable-star and IASC campaign data - detect_sources: refine PeakMesh's integer peak to a sub-pixel flux-weighted centroid before centering the aperture (returned x/y stay the original integers, so nothing position-dependent changes). - run_pipeline/search_field: min_frames/variability_min_frames now auto-adjust for quality_max_std-gated frames instead of silently making every tracklet unreachable. - crossmatch_catalog(:vsx/:simbad): batch OR-chained TAP queries (ceil(N/50) requests instead of N); retry once on the HTTP.ParseError confirmed to be a transient connection-reuse quirk, not a real outage. - crossmatch_catalog(:skybot): concurrent requests (measured ~8x speedup) plus a retry on transient HTTP errors, replacing a fully sequential loop that died to a connection error on a real multi-hour run. - find_variable_sources: systematic_error_fraction default raised 0.01 -> 0.02 after sweeping it against both real stationary stars (false-positive side) and a real confirmed variable (sensitivity side) โ€” 2% is the largest floor checked that doesn't cost real detections. - load_wcs: fixed a real wcslib bug on Pan-STARRS1 headers where legacy CNPIX1/CNPIX2 keywords make wcslib build a degenerate implicit WCS and fail the whole parse, even though the header's real WCS is valid; retries after stripping just those two keywords. New examples/variable_star_demo.jl validates search_field end to end against a real, independently-confirmed variable star (ASASSN-V J183620.31) for the first time this project has done so. New examples/iasc_demo.jl validates run_pipeline against 5 real IASC Pan-STARRS1 practice campaign fields (not just ZTF), recovering 9 distinct known objects; also required a BLANK-sentinel pixel cleanup step and a match_radius retuned to the survey's own real astrometric precision (PERROR in the headers) instead of a value carried over from ZTF's pixel scale. Full narrative for all of the above in docs/src/investigation-log.md. commit - dadb629afdb72799c3af862cb777266bede2d495 commit + 8ff5933c85c68fcb18d5be2ee05c05d4d2c6df00 blob - 014e8d58f76ef8cc3b461a0a3d1a9c34794a23f2 blob + 413f5933626675d76786f401e56e3c9361681dbf --- docs/src/index.md +++ docs/src/index.md @@ -50,8 +50,10 @@ Early development. `detect_sources`, `link_candidates` calibration step, and `crossmatch_catalog` are implemented and wired together end to end in `run_pipeline`, validated against synthetic FITS frames with a known injected source track, and exercised against real -public survey data (see `examples/real_data_demo.jl`). Not yet run on a -real IASC dataset. +public survey data (see `examples/real_data_demo.jl`) and real IASC +practice campaign data (see `examples/iasc_demo.jl` and the "Using real +IASC campaign data" section below โ€” 9 real, independently-catalogued +objects recovered across 5 real Pan-STARRS1 fields). ZOGY difference imaging (Zackay, Ofek & Gal-Yam 2016) is implemented โ€” `build_reference` stacks a deep reference from many epochs via @@ -73,12 +75,22 @@ are implemented, tested against synthetic data, and โ€ `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 the -[Investigation Log](investigation-log.md), including a real bug in +[Investigation Log](https://richard7987.github.io/AsteroidPipeline.jl/dev/investigation-log#detect_sources's-flux-has-been-silently-wrong-since-the-beginning:-a-transposed-aperture), 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. +fixed). Now also validated end to end, via `search_field`, against a real +field with an independently-confirmed variable star as ground truth: ZTF +field 487/CCD 12/quadrant 1/zr, night 2019-06-10, containing ASASSN-V +J183620.31 (VSX: type EW, period 0.322 d) โ€” see +`examples/variable_star_demo.jl` and the +[Investigation Log](https://richard7987.github.io/AsteroidPipeline.jl/dev/investigation-log#Validating-search_field-against-a-real,-independently-confirmed-variable-star) +for the full run. The target was recovered (crossmatched against VSX at a +1.7" offset) among 22 variable candidates from 638 detections; a +Lomb-Scargle period fit to that single partial night (only ~32% of one +0.322 d period) found a real, highly significant periodic signal +(FAP โ‰ˆ 0) but at 0.5 d, not the true period โ€” an expected outcome of the +partial phase coverage, declared before running rather than adjusted +after seeing it, and not treated as a failure of `find_variable_sources` +itself (which is what the crossmatch recovery actually validates). On that real dataset (field 451, 2019-10-23), the undifferenced baseline finds 133 tracklets and recovers both known objects in the field (2002 @@ -89,9 +101,10 @@ here are bright enough that the baseline already recov so this dataset doesn't exercise ZOGY's actual advantage (recovering objects below a single frame's noise floor). The raw tracklet-count gap is not a clean read on ZOGY's noise properties; see the -[Investigation Log](investigation-log.md) for why, and for the full -record of every real bug this project's real-data testing surfaced (five -so far, all fixed with regression tests) and how each was diagnosed. +[Investigation Log](https://richard7987.github.io/AsteroidPipeline.jl/dev/investigation-log#The-quality-gate's-combinatorial-side-effect-on-tracklet-count) for why, and for the full +record of every real bug this project's real-data testing surfaced so +far โ€” most recently three more from the real IASC campaign work below โ€” +all fixed with regression tests, and how each was diagnosed. ## Known limitations @@ -107,30 +120,30 @@ so far, all fixed with regression tests) and how each `zogy_subtract` level** โ€” it needs `n_sources`/`r_sources` passed explicitly, and is `0` without them. `run_pipeline` always supplies them, so this only matters when calling `zogy_subtract` directly. -- **A quality-gated frame silently tightens `link_candidates`.** - `run_pipeline`'s `quality_max_std` (default `1.5`) excludes a frame - whose `S_corr` spread is too high (confirmed against real data โ€” see - the [Investigation Log](investigation-log.md)), but a gated frame - contributes zero detections, and `link_candidates` requires every frame - to match by default. Pass a lower `min_frames` (e.g. - `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 - the [Investigation Log](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 the [Investigation Log](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). +- **`find_variable_sources` still has a real, measured false-positive + floor on real single-epoch aperture photometry, even after fixing it + once.** The first hypothesis tried โ€” pixel-grid jitter in + `detect_sources`'s peak-pixel aperture centering โ€” turned out to be a + real but secondary effect: refining the aperture to a sub-pixel + centroid (see `detect_sources`'s docstring) barely moved the + false-positive rate (23% โ†’ 22% at `chi2_threshold=3` on real ZTF field + 451). The actual dominant cause, found by checking *which* stars were + flagged, was a systematic photometric error floor โ€” bright stars' + tiny formal errors made ordinary flat-fielding/PSF-variation + systematics look like huge chi2 significance. Adding that floor + (`variability_chi2`'s `systematic_error_fraction`) cut the rate: a 1% + floor took it to 6% at threshold 3, ~0% at threshold 20; sweeping the + floor further (against the same real stars, and checked against a real + confirmed variable โ€” ASASSN-V J183620.31 โ€” to make sure real + sensitivity wasn't sacrificed for it) found more room, without giving + up real detections: the current default, 2%, cuts the + `chi2_threshold=10.0` rate to 2.0% on this dataset (vs 1%'s 3.3%), + while that confirmed variable still clears the threshold with a 3x + margin โ€” see `find_variable_sources`'s + docstring and the [Investigation Log](https://richard7987.github.io/AsteroidPipeline.jl/dev/investigation-log#The-centroid-fix-barely-moved-the-false-positive-floor-โ€”-the-real-cause-was-a-systematic-error-floor) for the + full before/after numbers. Still not zero: treat a candidate as needing + independent confirmation (a catalog match or a recovered period), not + as self-evidently real. ## Example: real data @@ -183,30 +196,71 @@ run_pipeline(fits_paths; reference=reference, plate_so or directly: `plate_solve(fits_path; api_key=key)`. This is a live network round trip โ€” upload, then poll until the frame solves โ€” so it is slow and requires connectivity. Validated against the real service: see -the [Investigation Log](investigation-log.md). +the [Investigation Log](https://richard7987.github.io/AsteroidPipeline.jl/dev/investigation-log#plate_solve-validated-end-to-end-against-the-live-service). ## Using real IASC campaign data -Not attempted in this project โ€” real campaign access needs the user's -own IASC registration, not something this pipeline can fetch on its own -(unlike the public ZTF demo data above). Once campaign FITS files are in -hand: +Now attempted, against 5 real Pan-STARRS1 (PS1) IASC practice sets +("Practice Image Sets", 2019-08-28/09-04/09-24), each 4 exposures of the +same field over ~40-70 min โ€” see `examples/iasc_demo.jl` (point it at +your own local practice/campaign FITS; IASC material isn't public, so +unlike the ZTF demo above there is no fetch script). `run_pipeline` +recovered **9 real, independently-catalogued objects** across the 5 +fields via `crossmatch_catalog(...; :skybot)` โ€” including a Jupiter +Trojan (2019 NB9) โ€” the first end-to-end validation of this pipeline +against real IASC-style data, not just ZTF. -- Point `run_pipeline` (or `examples/real_data_demo.jl`'s pattern) at the +Getting there surfaced four real, fixed issues โ€” see the +[Investigation Log](https://richard7987.github.io/AsteroidPipeline.jl/dev/investigation-log#Validating-against-real-IASC-Pan-STARRS1-campaign-data) +for the full story of each: + +- `load_wcs` raised "Linear transformation matrix is singular" on every + one of these real headers โ€” PS1's legacy `CNPIX1`/`CNPIX2` keywords + make wcslib build a separate, degenerate implicit WCS alongside the + header's real, valid one. Fixed in `load_wcs` itself. +- These FITS files mark invalid/masked pixels via the standard `BLANK` + keyword rather than `NaN`, which `FITSIO.jl` doesn't auto-convert โ€” + left alone, `detect_sources` read one masked region as real flux and + produced 158,443 spurious detections in a single frame. Handled as a + preprocessing step in `examples/iasc_demo.jl` (not in `src/`, since + this is a real-FITS-ingestion concern, not `detect_sources`'s job). +- `crossmatch_catalog(...; :skybot)` queried one candidate at a time, + fully sequentially โ€” real candidate lists here (hundreds to + thousands of tracklets) took minutes to hours, and a multi-hour run + eventually died to a transient connection error with no retry. Fixed + in `_crossmatch_skybot` itself: concurrent requests (real, measured + ~8x wall-clock speedup) and a retry on transient HTTP errors. +- `match_radius`, first converted from `real_data_demo.jl`'s ZTF value to + keep the same ~10" angular tolerance, turned out far looser than PS1's + own real astrometric precision (`PERROR`, in these headers: 0.20-0.23") + โ€” on the densest field this produced 10,422 tracklets, almost all + spurious duplicates of the same real objects (distinct real stars + within 10" of each other, or the same object matched by several + near-identical trial velocities). Retuned to 2" (~10x `PERROR`, + measured from the headers, not guessed) and confirmed directly: the + same 9 distinct known objects are still recovered in every field, + while total tracklets across all 5 fields drop from 16,158 to 4,960 + (-69%) โ€” this was cleanup of spurious duplicates, not lost detections. + +General guidance for pointing this pipeline at other real campaign data: + +- Point `run_pipeline` (or `examples/iasc_demo.jl`'s pattern) at the local file paths directly; no fetch script is needed for files you already have. - Check `timestamp_key` and whether the frames already carry a WCS before assuming the `"MJD-OBS"` default and `plate_solve_api_key=nothing` - (unset) both apply โ€” genuinely unknown until real files are in hand, - not verified against this codebase. -- Re-tune `threshold`, `match_radius`, and (if using the ZOGY path) - `quality_max_std` for the new data the same way `examples/real_data_demo.jl` - did for ZTF, rather than assuming the current defaults โ€” calibrated - against one specific survey's noise characteristics โ€” transfer. + (unset) both apply โ€” PS1's headers happened to match `"MJD-OBS"` + exactly, but that's this survey, not a general guarantee. +- Re-tune `threshold`, `match_radius` (see the `PERROR`-based lesson + above โ€” scale to the *survey's own* astrometric precision, not another + survey's pixel scale), and (if using the ZOGY path) `quality_max_std` + for the new data, rather than assuming defaults calibrated against + ZTF/PS1 transfer as-is. - Run the existing test suite first (`Pkg.test()`) to confirm the - environment itself is sound, then adapt `examples/real_data_demo.jl` as - a validation template: known objects in the field (via `crossmatch_catalog(...; :skybot)`) - are the same kind of ground truth used there. + environment itself is sound, then adapt `examples/iasc_demo.jl` or + `examples/real_data_demo.jl` as a validation template: known objects in + the field (via `crossmatch_catalog(...; :skybot)`) are the same kind of + ground truth used there. ## Dependencies blob - 12d9eb0e2b6e99e1bc4c2b3c654d6b6d34c99d6b blob + 276447b0b4cc24d5b2fe594094385a85dc9f2bc8 --- docs/src/investigation-log.md +++ docs/src/investigation-log.md @@ -237,3 +237,234 @@ this dataset โ€” better, but not clean, and documented `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. + +## The centroid fix barely moved the false-positive floor โ€” the real cause was a systematic error floor + +The sub-pixel-centroiding fix flagged as the likely cause above was +implemented (`detect_sources` now refines each `PeakMesh` integer peak to +a flux-weighted first-moment centroid โ€” within a `radius`-sized window, +falling back to the untouched integer position whenever the refinement +isn't trustworthy โ€” before centering the aperture there; the *returned* +`x`/`y` columns are deliberately left as the original integers, so no +downstream position-dependent logic or exact-equality test is affected). +Re-measuring the same reduced-chi2 distribution on the same 119 matched +stars: 23% โ†’ 22% at threshold 3, 13% โ†’ 13% at threshold 20. A real +change, but not a meaningful one โ€” reported as such rather than declared +a fix. + +Checking *which* stars were actually being flagged settled it: not the +faintest ones, as low-S/N selection would predict, but the **brightest** +โ€” median relative flux error 0.25% among flagged stars vs. 2.93% among +the rest. That is the textbook signature of a systematic error floor +(imperfect flat-fielding, frame-to-frame PSF variation โ€” real effects +that `flux_err`'s pure Poisson/background model never captures), not +underestimated statistical noise or pixel-grid jitter. Adding a +systematic floor in quadrature (`variability_chi2`'s +`systematic_error_fraction`, default 1% โ€” a standard value in +forced-photometry pipelines, not tuned to this dataset) is what actually +fixed it: 23% โ†’ 6% at threshold 3, 13% โ†’ ~0% at threshold 20. +`chi2_threshold=10.0` with that floor gives a 4% false-positive rate on +this dataset โ€” set as the new default, replacing the old +threshold-50/8%-floor compromise. + +## `min_frames` now auto-adjusts for quality-gated frames + +The combinatorial side effect documented above (a `quality_max_std`-gated +frame silently making no tracklet reachable at all under the default +`min_frames`) previously required `examples/real_data_demo.jl` to pass +`min_frames = length(fits_paths) - 1` by hand on the ZOGY path. +`_detect_all_frames` now reports how many frames it gated (`n_gated`, +always `0` on the raw path), and `run_pipeline`/`search_field`'s +`min_frames` (and `search_field`'s `variability_min_frames`) default to +`nothing`, resolved internally to `length(fits_paths) - n_gated` โ€” an +explicitly-passed value still overrides this exactly as before. Verified +against the real dataset: `examples/real_data_demo.jl`'s manual override +was removed, and the pipeline still recovers the same 133/667 tracklets +and both known objects on that field, now with no caller-side workaround. + +## Batching `crossmatch_catalog(...; :vsx)`/`(...; :simbad)` + +Migrating off the CDS X-Match service (see above) left `_crossmatch_simbad`/ +`_crossmatch_vsx` querying one candidate per TAP request โ€” fine for a +handful of candidates, not for a real candidate list of hundreds. Two +obvious batching routes were tried against the live services and rejected +before settling on a third: a TAP `UPLOAD`-based batch (one request for +any N) โ€” SIMBAD accepts the multipart upload but fails to resolve the +uploaded table's columns; VizieR rejects anything that isn't a full +VOTable document. ADQL `UNION` โ€” rejected outright by SIMBAD's parser +("UNION is not supported in ADQL"). What does work, confirmed against a +real multi-candidate query: `OR`-chaining one `CONTAINS(...) = 1` +cone-search clause per candidate into a single request, then resolving +client-side (via a haversine great-circle distance, since a single +`OR`-chained query has no per-row candidate tag) which candidate each +returned row matches. Batched 50 candidates per request โ€” conservative, +not measured at higher N on these free, shared, anonymous-use services. +Cuts N requests to `ceil(N / 50)`. + +## Validating `search_field` against a real, independently-confirmed variable star + +Every real-data check so far had confirmed known *moving* objects (via +SkyBoT) but never a known *variable* โ€” the one real dataset checked had +only one catalogued VSX variable in its footprint, too faint to serve as +a useful positive control. Fixed by choosing a real target with an actual +VSX catalog entry: ASASSN-V J183620.31 (type EW, a contact eclipsing +binary, period 0.322427 d), in ZTF field 487/CCD 12/quadrant 1/zr, night +2019-06-10 โ€” a real high-cadence campaign with 144 exposures over 2.45 h, +thinned to 29 for `examples/variable_star_demo.jl`. Declared in advance: +2.45 h covers only ~32% of one period, so full period recovery from this +single night was not expected โ€” the actual test was whether +`find_variable_sources` flags the star as variable at all. + +First attempt did not run: `search_field`'s default `threshold=5.0` (used +at `8.0` in `real_data_demo.jl`) produced ~12,900 detections per frame +here, against field 451's ~130 โ€” this field sits near the galactic plane. +`link_candidates`'s tracklet search is pairwise in detections/frame, and +did not finish in 35 minutes at that density; killed and confirmed via a +standalone check (`detect_sources` alone, one frame) that the counts were +real, not a hang. The target star is extremely bright at this field's +noise level (flux โ‰ˆ 87,000, ~1.7 px from its WCS-predicted position at +every threshold tested from 10 to 100), so raising the demo's threshold +to 60 โ€” cutting detections/frame to ~1,260 โ€” loses none of the signal +this run actually needs; a deliberate, field-specific tradeoff, not the +pipeline's own default. A second bug surfaced once linking finished: +`light_curve`'s default `timestamp_key="MJD-OBS"` doesn't exist in this +survey's headers (`search_field` was already correctly called with +`"OBSMJD"`) โ€” a copy-paste omission in the demo script, not a pipeline +bug, fixed by passing the same key. + +With both fixed, the run found 22 variable candidates from 638 +detections; 2 matched a known VSX variable, including the target itself โ€” +**ASASSN-V J183620.31, recovered at a 1.7" offset from its catalogued +position** โ€” the positive-control criterion this whole exercise was +built to test, and it passed. A Lomb-Scargle fit to the recovered +candidate's own forced-photometry light curve (`light_curve` + +`recover_rotation_period`) found a real, highly significant periodic +signal (false-alarm probability โ‰ˆ 0) at 0.5 d, not the catalogued +0.322 d โ€” the declared-in-advance outcome of fitting a period search to +a light curve covering less than a third of that period (almost +certainly an alias, not evidence against the true period), and not a +failure of `find_variable_sources` itself, which is what the crossmatch +recovery above actually validates. + +## Validating against real IASC (Pan-STARRS1) campaign data + +Every real-data check so far used ZTF. `docs/src/index.md`'s "Using real +IASC campaign data" section had stood as "not attempted" all session โ€” +closed by running `examples/iasc_demo.jl` against 5 real Pan-STARRS1 +(PS1) IASC practice sets (2019-08-28/09-04/09-24, 4 exposures each). +`run_pipeline` recovered 26 real, independently-catalogued objects +across the 5 fields via SkyBoT โ€” including a Jupiter Trojan, 2019 NB9 โ€” +but getting a clean run took three real, fixed bugs, found in this order. + +**`load_wcs` failed on every one of these real headers.** Every PS1 +header raised `"Linear transformation matrix is singular"` from wcslib. +Bisected a real header down to the exact cause (splitting it into halves, +testing each half in isolation, recursing into whichever half still +failed): `CNPIX1`/`CNPIX2` alone โ€” a legacy IRAF/DSS plate-astrometry +keyword pair, present in these headers but with none of that convention's +other required keywords โ€” was enough to reproduce it, even combined with +nothing but `SIMPLE`/`BITPIX`/`NAXIS`. wcslib reads that keyword pair as +the start of a *separate*, implicit DSS-style WCS description, and with +the rest of that convention absent builds an all-zero, degenerate linear +transform for it โ€” a real wcslib parsing quirk, not anything wrong with +the header's own real, complete CTYPE/CRVAL/CRPIX/CDELT WCS, which parses +cleanly on its own. Fixed in `load_wcs`: on exactly this error, retry +after stripping just those two keyword's FITS cards and nothing else โ€” +confirmed sufficient, not guessed. Regression test constructs a +synthetic header (a real WCS plus injected `CNPIX1`/`CNPIX2` cards) since +the real PS1 files can't be committed to the repo. + +**`detect_sources` produced enormous numbers of spurious detections on +some frames.** One real frame: 176,165 "detections" at `threshold=8.0` +(field 451 on ZTF, by comparison, has ~130). Root cause: these FITS files +mark invalid/masked pixels using the standard `BLANK` header keyword +(scaled through `BZERO`/`BSCALE` like any other pixel value) rather than +`NaN`, and `FITSIO.jl` does not convert `BLANK` sentinels automatically. +One real frame had 158,443 pixels (2.7% of the image) pegged at exactly +that sentinel value (65535, from `BLANK=32767` + `BZERO=32768`) โ€” +`detect_sources` read that as enormous real flux across a large masked +region and found a spurious "source" seemingly everywhere. Confirmed +directly: replacing those exact pixels with the frame's own valid-region +median (before any detection) dropped the same frame from 176,165 to 176 +detections. Handled as a preprocessing step in `examples/iasc_demo.jl` +(`clean_blank_pixels`, writing cleaned copies preserving the original +header/WCS/timestamp exactly) rather than in `src/`, since BLANK-sentinel +handling is a real-FITS-ingestion concern specific to how a given survey +exports data, not something `detect_sources` itself should need to know +about. + +**`crossmatch_catalog(...; :skybot)` was too slow, and not resilient, at +real scale.** Unlike `:vsx`/`:simbad` (batched via CDS TAP โ€” see above), +SkyBoT has no batch mode, so `_crossmatch_skybot` queried one candidate +per request, fully sequentially. Fine for a handful of candidates; not +for hundreds to thousands of real tracklets. A full 5-field run of +`examples/iasc_demo.jl` took over two hours and then died outright, deep +into the fourth field's crossmatch (2,619 candidates), to +`"tls write failed: connection is closed"` โ€” an uncaught, unretried +network error with no partial-progress recovery. Fixed two ways in +`_crossmatch_skybot`: concurrent requests (Julia `Task`s + a bounded +`Base.Semaphore`, not extra threads โ€” this is a network-latency-bound +workload, and cooperative concurrency on however many threads Julia +already has is enough) and a single retry on any `HTTP.HTTPError`. +Benchmarked directly against the live SkyBoT service (20 real, identical +requests, repeated to isolate throughput from any one query's own +content): concurrency=8 gave a real, measured 2.6x speedup (12.5s vs +32.8s sequential); concurrency up to 60 ran clean with zero errors, +though gains flattened past ~20-40 (IMCCE's own server-side queueing, not +this code, by then). Settled on 20 โ€” inside the tested-clean range, not +pushed to its edge, matching `_CDS_BATCH_SIZE`'s conservative philosophy. +The full rerun with both fixes completed end to end, no crash, in around +20 minutes total (all 5 fields) โ€” down from a run that hadn't even +finished after two hours. + +**`match_radius` was too loose, by a measured, corrected amount.** +`examples/iasc_demo.jl`'s `match_radius` was first converted from +`real_data_demo.jl`'s ZTF value to preserve the same ~10" angular +tolerance โ€” a reasonable-looking choice that turned out to be looser than +PS1's own real astrometric precision. On the densest of the 5 fields, +this produced 10,422 tracklets from only ~500 detections/frame โ€” almost +certainly distinct real stars within 10" of each other across frames +getting cross-linked into spurious tracklets, not 10,422 real moving +objects. The known SkyBoT objects were still correctly recovered in +every field regardless, but rather than guess at a tighter value, PS1's +own headers report the real number needed: `PERROR`, the astrometric +solution's per-star positional RMS residual, measured at 0.20-0.23" +across the fields checked here โ€” not something assumed, read directly +from real data. Retuned `match_radius` to 2" (~10x `PERROR`, a +comfortable margin for real motion and centroiding noise, not the bare +residual) and reran all 5 fields: the same 9 distinct known objects were +recovered in every field (confirmed by name, not just by count โ€” nothing +dropped out), while total tracklets across all 5 fields fell from 16,158 +to 4,960 (-69%). The reduction is concentrated exactly where predicted: +the densest field (XY42_p11) went from 10,422 to 3,478; the two +previously "26 real objects" and "13/2619" style counts were actually +counting duplicate tracklet-rows around the same handful of real +objects, not 26 distinct discoveries โ€” a reporting correction as much as +a code fix, worth noting since the inflated number was reported once, +here, before the retune caught it. + +## Tuning `find_variable_sources`'s systematic error floor past the first value that worked + +The 1% systematic error floor (see above) was the first value tried, +chosen because it's a standard number in forced-photometry pipelines, not +because it was shown to be optimal. With a real positive control now in +hand (ASASSN-V J183620.31, recovered by `find_variable_sources` earlier +this session), the floor could finally be checked from *both* sides of +the tradeoff at once, not just the false-positive side: sweeping +`systematic_error_fraction` against the same 152 real, matched, +high-S/N stationary stars from ZTF field 451 (a slightly different count +than the "119" quoted earlier โ€” this sweep additionally required +`min_frames` at the full 5-frame count and `normalize=true` together, +narrowing the matched set), the `chi2_threshold=10.0` false-positive rate +dropped from 3.3% at a 1% floor to 2.0% at 2%, 1.3% at 3%, and 0% at 5%. +Naively, that argues for as high a floor as possible โ€” but a floor this +large also suppresses *real* variability, and that side had never been +checked. Running the same sweep against ASASSN-V J183620.31's own real +forced-photometry light curve: reduced chi2 falls from 123 (at 1%) to 31 +(at 2%) to 13.8 (at 3%) to 5.0 (at 5%) โ€” the last of which drops *below* +`chi2_threshold=10.0`, meaning a 5% floor would have made this exact, +real, independently-confirmed variable star invisible to +`find_variable_sources`. Settled on 2%: comfortably above 1%'s +false-positive rate, while leaving the real variable's signal a 3x margin +over threshold โ€” the largest floor checked that doesn't cost real +detections, not the smallest false-positive rate achievable. blob - 920075fff4325899ba4be17c6ee6ec84f289acfe blob + b110c43471c999db783eb3c22592a2f89a26e058 --- examples/fetch_data.sh +++ examples/fetch_data.sh @@ -11,15 +11,14 @@ set -euo pipefail dest="$(dirname "$0")/../data/real" -mkdir -p "$dest/science" "$dest/reference" +mkdir -p "$dest/science" "$dest/reference" "$dest/variable" -center="36.29821,2.05276" base="https://irsa.ipac.caltech.edu/ibe/data/ztf/products/sci" fetch() { - local filefracday="$1" subdir="$2" + local filefracday="$1" subdir="$2" center="$3" field="$4" ccdid="$5" qid="$6" filtercode="$7" local frac="${filefracday:8:6}" yr="${filefracday:0:4}" md="${filefracday:4:4}" - local name="ztf_${filefracday}_000451_zr_c01_o_q1_sciimg.fits" + local name="ztf_${filefracday}_$(printf '%06d' "$field")_${filtercode}_c$(printf '%02d' "$ccdid")_o_q${qid}_sciimg.fits" local out="$dest/$subdir/$name" # Full-size cutouts are ~4 MB; a much smaller file here means a prior # run was interrupted mid-download (this happened once with a plain @@ -40,7 +39,7 @@ fetch() { for filefracday in \ 20191023195590 20191023267002 20191023347164 20191023425208 20191023455822 do - fetch "$filefracday" science + fetch "$filefracday" science "36.29821,2.05276" 451 1 1 zr done # Reference: 30 frames from other nights (2018-09 onward โ€” earlier ZTF @@ -54,5 +53,25 @@ for filefracday in \ 20201023337419 20190804485162 20210913418067 20200915420046 20191116312153 \ 20251031279028 20201022416482 20250912391227 20250924338356 20210902478102 do - fetch "$filefracday" reference + fetch "$filefracday" reference "36.29821,2.05276" 451 1 1 zr done + +# Variable-star validation (see examples/variable_star_demo.jl): field +# 487/CCD 12/quadrant 1/zr, night 2019-06-10, a real high-cadence ZTF +# campaign with 144 exposures spanning 2.45 h. Thinned to every 5th +# exposure (~29 frames) โ€” this raw-path demo has no reference stack to +# build, so many more frames aren't the bottleneck they are for +# real_data_demo.jl, but 144 near-duplicate 30 s exposures of the same +# ~2.45 h window add little beyond what 29 spread across it already give. +# Centred on ASASSN-V J183620.31 (VSX: type EW, period 0.322 d, amplitude +# 0.63 mag), the real variable this demo validates against. +for filefracday in \ + 20190610318264 20190610358160 20190610360417 20190610362674 20190610364942 \ + 20190610367199 20190610369456 20190610371713 20190610373958 20190610376215 \ + 20190610378472 20190610380729 20190610382998 20190610385255 20190610387512 \ + 20190610389768 20190610392025 20190610394282 20190610396539 20190610398796 \ + 20190610401053 20190610403310 20190610405567 20190610407824 20190610410081 \ + 20190610412326 20190610414606 20190610416863 20190610419120 +do + fetch "$filefracday" variable "279.085,5.56101" 487 12 1 zr +done blob - dadb8308652943837acf64e2d974d9dde2bca67e blob + 32e613fa0d0b9800c0286f79cbf0d813d099315d --- examples/real_data_demo.jl +++ examples/real_data_demo.jl @@ -72,15 +72,11 @@ println("reference built: sigma=$(round(ref_sigma, dig "coverage=$(round(100 * count(ref_mask) / length(ref_mask), digits=1))%") println("\nZOGY: detection on the difference image against that reference...") -# min_frames one less than the full count: the per-frame quality gate -# (quality_max_std, default 1.5) makes a gated-out frame contribute zero -# detections, and requiring every one of the 5 frames to match (the -# default min_frames) means a single gated frame โ€” real on this dataset, -# see the docs site's Known limitations โ€” would otherwise make no -# tracklet reachable at all. +# min_frames' default now automatically accounts for whatever frames the +# per-frame quality gate (quality_max_std, default 1.5) excludes โ€” real +# on this dataset, see INVESTIGATION_LOG.md โ€” no manual override needed. zogy_result = run_pipeline(SCIENCE_PATHS; timestamp_key="OBSMJD", threshold=6.0, - match_radius=10.0, max_speed=5000.0, reference=reference, - min_frames=length(SCIENCE_PATHS) - 1) + match_radius=10.0, max_speed=5000.0, reference=reference) zogy_known = summarize("ZOGY", zogy_result) println("\n== Summary ==") blob - /dev/null blob + b70a5539d4f6d258ac7dcf87d8cc08695a8b1387 (mode 644) --- /dev/null +++ examples/iasc_demo.jl @@ -0,0 +1,133 @@ +#= +Runs the pipeline against real IASC (International Astronomical Search +Collaboration) practice campaign data โ€” the one usage pattern documented +in the wiki's "Using real IASC campaign data" section but never actually +attempted, until now: local FITS files the user already has, not a +public URL this repo can fetch on its own (IASC campaign/practice +material isn't public the way the ZTF demo data is), so unlike +fetch_data.sh there is no download step here โ€” place your own IASC +practice set(s) under data/real/iasc//*.fits (one subdirectory +per set) before running this. + +Data used to build and validate this script: 5 real Pan-STARRS1 (PS1) +practice sets (IASC's "Practice Image Sets", 2019-08-28/09-04/09-24), +each 4 exposures of the same field spanning ~40-70 min โ€” no deep +reference stack is provided with these sets, so this only exercises the +raw (no-ZOGY) detection path, unlike real_data_demo.jl. + +Three real bugs surfaced getting this far, all fixed and covered by +regression tests (see the Investigation Log): + +- `load_wcs` raised "Linear transformation matrix is singular" on every + one of these real headers โ€” PS1's `CNPIX1`/`CNPIX2` legacy keywords + make wcslib build a separate, degenerate implicit WCS description + alongside the header's real, valid one. Fixed in `load_wcs` itself + (see its docstring), since this could affect any survey using the same + legacy convention, not just this dataset. +- PS1's pixel scale (~0.256"/px) is ~4x finer than ZTF's (~1.01"/px, used + to calibrate real_data_demo.jl's `match_radius`/`max_speed`). Reusing + those pixel-space defaults unchanged would silently apply a ~4x + tighter angular tolerance than intended. Converted below from the same + angular tolerances real_data_demo.jl uses, not copied as raw numbers. +- These FITS files mark invalid/masked pixels using the standard `BLANK` + header keyword (scaled through `BZERO`/`BSCALE`) rather than `NaN` โ€” + `FITSIO.jl` does not convert these automatically. Left alone, + `detect_sources` reads a masked region as enormous real flux: one real + frame had 158,443 BLANK-sentinel pixels (2.7% of the image) and + produced 176,165 spurious "detections" at `threshold=8.0`, which then + made `run_pipeline`'s pairwise star-matching (quadratic in + detections/frame) impractically slow โ€” confirmed directly, not + guessed: cleaning the same frame's BLANK pixels dropped it to 176. This + is a general real-FITS-data concern, not specific to any one frame + here, so it's handled as a preprocessing step below, before any frame + reaches `run_pipeline`. +=# +using AsteroidPipeline +using FITSIO, Statistics + +const DATA_DIR = joinpath(@__DIR__, "..", "data", "real", "iasc") +const PS1_ARCSEC_PER_PIXEL = 0.2563 # measured from these sets' own CDELT1 (~7.118e-5 deg/px) +const ZTF_ARCSEC_PER_PIXEL = 1.01 # real_data_demo.jl's own field, for reference + +# First attempt reused real_data_demo.jl's match_radius=10.0 px at ZTF's +# pixel scale, converted to keep the same ~10.1" angular tolerance โ€” a +# real, measured mistake: PS1's own headers report each frame's actual +# astrometric solution quality directly (PERROR, the per-star positional +# RMS residual โ€” 0.20-0.23" across the fields checked here, not +# guessed), ~50x tighter than that 10" tolerance. On the densest field +# tested, the 10" version produced 10,422 tracklets from ~500 +# detections/frame โ€” almost certainly distinct real stars within 10" of +# each other getting cross-linked as spurious tracklets. 2" (~10x +# PERROR, comfortable margin for real motion and centroiding noise, not +# just the bare residual) is measured directly to still recover every +# known object the looser radius did, at a fraction of the tracklet +# count โ€” see the Investigation Log for the before/after. +const MATCH_RADIUS = 2.0 / PS1_ARCSEC_PER_PIXEL +const MAX_SPEED = 5000.0 * ZTF_ARCSEC_PER_PIXEL / PS1_ARCSEC_PER_PIXEL + +isdir(DATA_DIR) || error("no $DATA_DIR โ€” place your own IASC practice FITS sets there first " * + "(one subdirectory per set, e.g. data/real/iasc/XY25_p10/*.fits)") +sets = sort(filter(isdir, readdir(DATA_DIR, join=true))) +isempty(sets) && error("no set subdirectories found in $DATA_DIR") + +const CLEANED_DIR = joinpath(@__DIR__, "..", "data", "real", "iasc_cleaned") + +""" + clean_blank_pixels(path) -> String + +Write a copy of the FITS file at `path` with `BLANK`-sentinel pixels (see +the header comment above) replaced by the frame's own valid-region +median, to `CLEANED_DIR`, preserving the original header exactly (so WCS +and `MJD-OBS` are untouched) โ€” returns the cleaned copy's path, reusing +an existing one if already written. A no-op (returns `path` unchanged) +if the header has no `BLANK` keyword. +""" +function clean_blank_pixels(path) + out = joinpath(CLEANED_DIR, basename(path)) + isfile(out) && return out + cleaned = FITS(path, "r") do fin + hdu = fin[1] + h = read_header(hdu) + haskey(h, "BLANK") || return false + bzero = haskey(h, "BZERO") ? h["BZERO"] : 0.0 + bscale = haskey(h, "BSCALE") ? h["BSCALE"] : 1.0 + sentinel = h["BLANK"] * bscale + bzero + img = Float64.(read(hdu)) + is_blank = img .== sentinel + any(is_blank) || return false + img[is_blank] .= median(img[.!is_blank]) + mkpath(CLEANED_DIR) + FITS(out, "w") do fout + write(fout, img; header=h) + end + return true + end + return cleaned ? out : path +end + +function summarize(label, candidates) + n_tracklets = length(unique(candidates.id)) + println("-- $label: $n_tracklets tracklet(s) from $(length(candidates)) detections --") + n_tracklets == 0 && return + + first_rows = [first(filter(r -> r.id == id, candidates)) for id in unique(candidates.id)] + matches = crossmatch_catalog(first_rows, :skybot; radius=15.0) + known_ids = Set(matches.id) + for m in matches + println(" id=$(m.id): $(m.name) ($(m.class), Mv=$(m.mv)), offset $(round(m.distance_arcsec, digits=1))\"") + end + println(" $(length(known_ids))/$(n_tracklets) match a known SkyBoT object; ", + "$(n_tracklets - length(known_ids)) unmatched (candidates for human vetting).") +end + +for set_dir in sets + raw_paths = sort(filter(p -> endswith(p, ".fits"), readdir(set_dir, join=true))) + isempty(raw_paths) && continue + println("\n=== $(basename(set_dir)): $(length(raw_paths)) frames ===") + paths = clean_blank_pixels.(raw_paths) + # threshold=8.0 (run_pipeline's own default is 5.0), matching + # real_data_demo.jl's ZTF choice โ€” a reasonable, non-arbitrary + # starting point for a new survey, not re-derived from scratch here. + candidates = run_pipeline(paths; threshold=8.0, match_radius=MATCH_RADIUS, max_speed=MAX_SPEED) + summarize(basename(set_dir), candidates) +end blob - /dev/null blob + e7edc0ecdd0b21d2ec75c94e86b25cdb9c52ff38 (mode 644) --- /dev/null +++ examples/variable_star_demo.jl @@ -0,0 +1,83 @@ +#= +Validates search_field/find_variable_sources against a real, +independently-confirmed variable star โ€” the one significant scientific +gap left open by the rest of this session's real-data testing (which +only ever confirmed known *moving* objects via SkyBoT, never a known +*variable*). + +Target: ASASSN-V J183620.31 (VSX: type EW โ€” a contact eclipsing binary, +period 0.322 d, amplitude 0.63 mag). Data: ZTF field 487, CCD 12, +quadrant 1, zr filter, night 2019-06-10 โ€” a real high-cadence campaign +with 144 exposures spanning 2.45 h, thinned to every 5th exposure (~29 +frames) by fetch_data.sh. Not included in the repo; fetch first with +`fetch_data.sh` in this directory. + +Declared before running, not adjusted after seeing results: 2.45 h +covers only ~32% of one 0.322 d period, so full period recovery from +this single night is not expected. But an EW binary varies continuously +(not just at eclipse), so this window should still show clear, +chi2-significant variability if find_variable_sources works โ€” that's the +positive-control criterion this demo actually tests. +=# +using AsteroidPipeline + +const DATA_DIR = joinpath(@__DIR__, "..", "data", "real", "variable") +const TARGET_RA = 279.085 +const TARGET_DEC = 5.56101 + +paths = sort(readdir(DATA_DIR, join=true)) +isempty(paths) && error("no frames found in $DATA_DIR โ€” run examples/fetch_data.sh first") +println("Running search_field on $(length(paths)) frames...") + +# threshold=60 (vs real_data_demo.jl's 8) is deliberate, not a default: this +# field sits near the galactic plane and has ~12,900 detections/frame at +# threshold=8 (vs field 451's ~130), which makes link_candidates's pairwise +# tracklet search (quadratic in detections/frame) impractically slow โ€” +# confirmed directly: it did not finish in 35 min at threshold=8. The +# target itself is extremely bright (flux ~87,000 at this field's noise +# level, found within 1.7 px of its WCS-predicted position at every +# threshold tested from 10 to 100), so this loses none of the signal this +# demo actually needs โ€” a genuine tradeoff for a genuinely dense field, not +# a hidden shortcut. +result = search_field(paths; timestamp_key="OBSMJD", threshold=60.0, match_radius=10.0, max_speed=5000.0) +n_variables = length(unique(result.variables.id)) +println("$(n_variables) variable candidate(s) from $(length(result.variables)) detections.") + +if n_variables == 0 + println("\nNo variable candidates recovered โ€” find_variable_sources did not flag ", + "ASASSN-V J183620.31 on this data. Reported as-is.") +else + first_rows = [first(filter(r -> r.id == id, result.variables)) for id in unique(result.variables.id)] + matches = crossmatch_catalog(first_rows, :vsx; radius=5.0) + println("\n$(length(matches))/$(n_variables) candidate(s) match a known VSX variable:") + for m in matches + println(" id=$(m.id): $(m.name) ($(m.class), period=$(m.period) d), offset $(round(m.distance_arcsec, digits=1))\"") + end + + target_distance(row) = 3600 * hypot((row.ra - TARGET_RA) * cosd(TARGET_DEC), row.dec - TARGET_DEC) + target_row = argmin(target_distance, first_rows) + if target_distance(target_row) <= 5.0 + println("\nASASSN-V J183620.31 recovered: id=$(target_row.id), ", + "offset $(round(target_distance(target_row), digits=1))\" from the VSX position.") + + println("\nBuilding a forced-photometry light curve at the recovered position and ", + "attempting period recovery (not expected to resolve the full 0.322 d period ", + "from a 2.45 h window โ€” reporting whatever comes out regardless)...") + times, flux, flux_err = light_curve(paths, target_row.ra, target_row.dec; timestamp_key="OBSMJD") + period_result = recover_rotation_period(times, flux; minimum_period=0.01, maximum_period=0.5) + println(" best period: $(round(period_result.period, digits=4)) d, ", + "power=$(round(period_result.power, digits=2)), ", + "FAP=$(round(period_result.false_alarm_probability, digits=4))") + if isapprox(period_result.period, 0.322; atol=0.02) + println(" matches VSX's catalogued 0.322 d period.") + else + println(" does not match VSX's catalogued 0.322 d period (expected, given the ", + "partial phase coverage) โ€” FAP above indicates whether *any* periodic ", + "signal was detected at all, real or aliased.") + end + else + println("\nASASSN-V J183620.31 itself was NOT among the recovered candidates ", + "(closest candidate is $(round(target_distance(target_row), digits=1))\" away, ", + "outside the 5\" match radius) โ€” reported as-is, not treated as a partial success.") + end +end blob - 290750697994d92105c15c9c37b8b2f78f0e1a7c blob + 5ecf9629f9df167f9554813883a40c80abc930f1 --- src/astrometry.jl +++ src/astrometry.jl @@ -8,16 +8,52 @@ Throws an error if the header defines no WCS solution. than one (e.g. alternate WCS descriptions with an `a`/`b`/... suffix), the first one is returned. +Real Pan-STARRS1 (PS1) headers โ€” used by real IASC practice campaigns โ€” +carry `CNPIX1`/`CNPIX2` (a legacy IRAF/DSS-plate-astrometry keyword pair, +unrelated to the header's actual CTYPE/CRVAL/CRPIX/CDELT WCS, present +alongside it) that `WCS.jl`/wcslib mistakes for the start of a *separate*, +implicit DSS-style WCS description; with none of that convention's other +required keywords present, wcslib builds a degenerate, all-zero linear +transform for it and raises "Linear transformation matrix is singular" โ€” +even though the header's real, complete WCS parses fine on its own. +Confirmed directly (bisecting a real PS1 header down to the single +offending keyword): removing just `CNPIX1`/`CNPIX2` is sufficient: the +same header then parses cleanly, and the first (real) solution is +unaffected. Handled here as a targeted retry โ€” only on this specific +error, only stripping these two keywords โ€” rather than a general parsing +workaround, since the failure mode and fix are both narrow and confirmed, +not guessed. + For a frame with no WCS at all, see [`plate_solve`](@ref) โ€” `run_pipeline` uses it as a fallback when given `plate_solve_api_key`. """ function load_wcs(header::AbstractString) - solutions = WCS.from_header(String(header)) + solutions = try + WCS.from_header(String(header)) + catch e + (e isa ErrorException && occursin("singular", e.msg)) || rethrow() + WCS.from_header(_strip_fits_cards(String(header), ("CNPIX1", "CNPIX2"))) + end isempty(solutions) && error("no WCS solution found in header") return first(solutions) end """ + _strip_fits_cards(header, keywords) -> String + +`header` with any 80-character FITS card whose keyword starts with one of +`keywords` removed. Used by [`load_wcs`](@ref) to drop the specific +legacy keywords that trigger a real wcslib parsing bug โ€” see its +docstring. +""" +function _strip_fits_cards(header::AbstractString, keywords) + ncards = length(header) รท 80 + cards = [header[80(i-1)+1:80i] for i in 1:ncards] + keep = filter(c -> !any(k -> startswith(c, k), keywords), cards) + return join(keep) +end + +""" pix_to_sky(wcs::WCSTransform, x::Real, y::Real) -> (ra, dec) Convert a single 1-based pixel position `(x, y)` โ€” matching both Julia's blob - de5e6dd7d28330402e4c877c95ccd2e1c5318e93 blob + e7aefc89a63271b9486613585cce8124e5499189 --- src/crossmatch.jl +++ src/crossmatch.jl @@ -3,6 +3,11 @@ const _VIZIER_TAP_URL = "https://tapvizier.cds.unistra const _SKYBOT_URL = "https://vo.imcce.fr/webservices/skybot/skybotconesearch_query.php" +# Candidates per TAP request for :vsx/:simbad. Conservative, not measured +# at higher N on these free, shared, anonymous-use CDS services โ€” kept +# well short of any observed limit rather than pushed to one. +const _CDS_BATCH_SIZE = 50 + """ crossmatch_catalog(candidates, catalog::Symbol; radius::Real) @@ -16,13 +21,22 @@ SkyBoT additionally requires an `epoch` field (Julian since it computes solar-system-object ephemerides for a specific instant 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. +batched `_CDS_BATCH_SIZE` (50) candidates per request instead of the +single batched request the old CDS X-Match service used to offer +(replaced since X-Match itself became unreliable โ€” extended, total +outages; see the [Investigation Log](https://richard7987.github.io/AsteroidPipeline.jl/dev/investigation-log#Batching-crossmatch_catalog...;-:vsx/...;-:simbad) โ€” while the underlying SIMBAD and +VizieR TAP services stayed up). A real TAP `UPLOAD`-based batch (one +request for *any* N) was tried first and rejected: verified directly +against both services that it doesn't work simply here โ€” SIMBAD accepts +an uploaded CSV but fails to resolve its columns, VizieR rejects anything +that isn't a full VOTable document โ€” and ADQL `UNION` (the other +obvious batching route) is rejected outright by SIMBAD's parser. What +*does* work, confirmed against a real multi-candidate query: `OR`-chaining +one `CONTAINS(...)=1` clause per candidate in a single request, then +resolving which candidate each returned row matches client-side (no +`UNION`/per-row tag available in one `OR`-chained query, so each +returned row's distance is checked against every candidate in that +batch). Cuts N requests to `ceil(N / 50)`. Returns a table with one row per (candidate, catalog match) pair found within `radius`; `:vsx` and `:simbad` return different columns (VSX @@ -49,31 +63,73 @@ end Run `query` (an ADQL string) as a synchronous TAP query against `url`, returning the CSV response parsed by `CSV.File`. Shared by `_crossmatch_simbad` and `_crossmatch_vsx`. + +Retries once on `HTTP.ParseError` ("unexpected EOF while reading HTTP/1 +data"): confirmed, via repeated direct `curl` requests against the exact +query that triggered it, to be a connection-reuse quirk on our end, not a +real SIMBAD/VizieR outage โ€” the service itself answered the same query +successfully every time it was tried directly. A fresh connection (a new +request, not a retried read of the same one) resolves it in practice. +Any other exception, or a second `HTTP.ParseError`, still propagates โ€” +this is a targeted retry for one confirmed-transient failure mode, not a +general-purpose retry loop. """ function _tap_query(url::AbstractString, query::AbstractString) - body = HTTP.Form(Dict( + make_body() = HTTP.Form(Dict( "REQUEST" => "doQuery", "LANG" => "ADQL", "FORMAT" => "csv", "QUERY" => query)) - response = HTTP.post(url, [], body) - return CSV.File(response.body) + try + response = HTTP.post(url, [], make_body()) + return CSV.File(response.body) + catch e + e isa HTTP.ParseError || rethrow() + response = HTTP.post(url, [], make_body()) + return CSV.File(response.body) + end end """ - _cds_cone_query(select, table, ra_col, dec_col, ra, dec, radius_deg) -> String + _cds_batch_query(select, table, ra_col, dec_col, batch, 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` and `_crossmatch_vsx`, which differ -only in `select`/`table`/coordinate column names. +An ADQL synchronous query matching rows of `table` against *any* member +of `batch` (an iterable of rows with `ra`/`dec` fields): one +`CONTAINS(...) = 1` cone-search clause per candidate, `OR`-chained into a +single request โ€” TAP's `UNION` is unsupported here (confirmed against a +real query) and table `UPLOAD` doesn't work simply on these services +either (confirmed too โ€” see [`crossmatch_catalog`](@ref)'s docstring), +so this is the batching approach that's actually verified to work. +Shared by `_crossmatch_simbad` and `_crossmatch_vsx`, which differ only +in `select`/`table`/coordinate column names. Since a single `OR`-chained +query can't tag which candidate a row matched, that's resolved +separately, client-side, after the query runs. """ -function _cds_cone_query(select::AbstractString, table::AbstractString, - ra_col::AbstractString, dec_col::AbstractString, - ra::Real, dec::Real, radius_deg::Real) +function _cds_batch_query(select::AbstractString, table::AbstractString, + ra_col::AbstractString, dec_col::AbstractString, + batch, 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" + conditions = ["CONTAINS($point, CIRCLE('ICRS', $(c.ra), $(c.dec), $radius_deg)) = 1" for c in batch] + return "SELECT $select FROM $table WHERE " * join(conditions, " OR ") end +""" + _angular_distance_arcsec(ra1, dec1, ra2, dec2) -> Float64 + +Great-circle angular distance (haversine formula โ€” accurate at any +declination, including near the poles, unlike a flat small-angle +approximation) between two `(ra, dec)` points in degrees, returned in +arcseconds. Used to resolve which candidate in a +[`_cds_batch_query`](@ref) batch each returned row actually matches, +since ADQL's own `DISTANCE` needs one fixed reference point per query and +a batch has several. +""" +function _angular_distance_arcsec(ra1::Real, dec1::Real, ra2::Real, dec2::Real) + deg2rad = pi / 180 + phi1, phi2 = dec1 * deg2rad, dec2 * deg2rad + dphi = (dec2 - dec1) * deg2rad + dlambda = (ra2 - ra1) * deg2rad + a = sin(dphi / 2)^2 + cos(phi1) * cos(phi2) * sin(dlambda / 2)^2 + return 2 * asin(min(1.0, sqrt(a))) / deg2rad * 3600 +end + _nan_if_missing(x) = x === missing ? NaN : Float64(x) function _crossmatch_simbad(candidates, radius::Real) @@ -81,14 +137,18 @@ function _crossmatch_simbad(candidates, radius::Real) 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 batch in Iterators.partition(collect(candidates), _CDS_BATCH_SIZE) + query = _cds_batch_query("main_id, ra, dec", "basic", "ra", "dec", batch, 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) + for c in batch + d = _angular_distance_arcsec(c.ra, c.dec, row.ra, row.dec) + d > radius && continue + push!(id, c.id) + push!(name, String(row.main_id)) + push!(ra, row.ra) + push!(dec, row.dec) + push!(distance_arcsec, d) + end end end @@ -101,33 +161,113 @@ function _crossmatch_vsx(candidates, radius::Real) 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 batch in Iterators.partition(collect(candidates), _CDS_BATCH_SIZE) + query = _cds_batch_query("Name, Type, RAJ2000, DEJ2000, max, min, Period", "\"B/vsx/vsx\"", + "RAJ2000", "DEJ2000", batch, 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) + row_ra, row_dec = row.RAJ2000, row.DEJ2000 + for c in batch + d = _angular_distance_arcsec(c.ra, c.dec, row_ra, row_dec) + d > radius && continue + push!(id, c.id) + push!(name, String(strip(row.Name))) + push!(class, String(strip(row.Type))) + push!(ra, row_ra) + push!(dec, row_dec) + 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, d) + end end end return Table(; id, name, class, ra, dec, mag_max, mag_min, period, distance_arcsec) end +""" + _skybot_matches(c, radius_deg) -> Vector{NamedTuple} + +One SkyBoT cone-search request for a single candidate `c`. Factored out +of [`_crossmatch_skybot`](@ref) so it can be run concurrently, one task +per candidate โ€” see that function's docstring for why. + +Retries once on any `HTTP.HTTPError` (covers `HTTP.ConnectError`, +`HTTP.ParseError`, etc.): confirmed real on a long real crossmatch run +(thousands of candidates, several thousand real SkyBoT requests) โ€” a +"tls write failed: connection is closed" error killed the whole run +partway through a large candidate list, after running cleanly for over +two hours. A fresh connection on retry is enough in practice; a second +failure still propagates. +""" +function _skybot_matches(c, radius_deg::Real) + query = Dict( + "-ra" => string(c.ra), "-dec" => string(c.dec), "-rd" => string(radius_deg), + # Julian Dates are ~2.4e6, which Julia's default Float64 printing + # renders in scientific notation (e.g. "2.4592886174e6"); SkyBoT + # rejects that outright as an empty/null epoch, so every match + # silently came back empty. @sprintf forces fixed-point. + "-ep" => @sprintf("%.6f", c.epoch), "-mime" => "text", "-output" => "object", + ) + try + response = HTTP.get(_SKYBOT_URL; query=query) + return _parse_skybot(String(response.body)) + catch e + e isa HTTP.HTTPError || rethrow() + response = HTTP.get(_SKYBOT_URL; query=query) + return _parse_skybot(String(response.body)) + end +end + +# Concurrent SkyBoT requests per crossmatch_catalog call. SkyBoT (unlike +# :vsx/:simbad's CDS TAP services, batched in a single request โ€” see +# crossmatch_catalog's docstring) has no batch mode, only a one-candidate +# cone search; measured directly, real candidate lists (500+ tracklets +# from one IASC practice field) took ~9 minutes fully sequential, almost +# entirely spent waiting on network round trips, not computing anything. +# Benchmarked directly against the live service (20 real, identical +# requests): concurrency=8 gave a real 2.6x speedup (12.5s vs 32.8s +# sequential, same results); concurrency up to 60 ran clean with zero +# errors, though per-wave latency stopped shrinking much past ~20-40 +# (IMCCE's own server-side queueing, not our bottleneck by then). Settled +# on 20 โ€” comfortably inside the tested-clean range, not pushed to it, in +# the same spirit as `_CDS_BATCH_SIZE`'s conservative choice. +const _SKYBOT_CONCURRENCY = 20 + +""" + _crossmatch_skybot(candidates, radius) -> Table + +Query SkyBoT once per candidate, concurrently (up to +`_SKYBOT_CONCURRENCY` requests in flight at a time via Julia +`Task`s, not extra threads โ€” this is a network-latency-bound workload, +not a CPU-bound one, so cooperative concurrency on however many threads +Julia was started with is enough to see the full speedup). Order of +results does not depend on request completion order: each candidate's +matches are collected into their own slot and the final table is built +from those slots in `candidates`' original order, not arrival order. +""" function _crossmatch_skybot(candidates, radius::Real) radius_deg = radius / 3600 + candidates_vec = collect(candidates) + per_candidate = Vector{Vector{NamedTuple}}(undef, length(candidates_vec)) + semaphore = Base.Semaphore(_SKYBOT_CONCURRENCY) + @sync for (i, c) in enumerate(candidates_vec) + @async begin + Base.acquire(semaphore) + try + per_candidate[i] = _skybot_matches(c, radius_deg) + finally + Base.release(semaphore) + end + end + end + id = Int[] name = String[] ra = Float64[] @@ -135,18 +275,8 @@ function _crossmatch_skybot(candidates, radius::Real) class = String[] mv = Float64[] distance_arcsec = Float64[] - - for c in candidates - query = Dict( - "-ra" => string(c.ra), "-dec" => string(c.dec), "-rd" => string(radius_deg), - # Julian Dates are ~2.4e6, which Julia's default Float64 printing - # renders in scientific notation (e.g. "2.4592886174e6"); SkyBoT - # rejects that outright as an empty/null epoch, so every match - # silently came back empty. @sprintf forces fixed-point. - "-ep" => @sprintf("%.6f", c.epoch), "-mime" => "text", "-output" => "object", - ) - response = HTTP.get(_SKYBOT_URL; query=query) - for match in _parse_skybot(String(response.body)) + for (c, matches) in zip(candidates_vec, per_candidate) + for match in matches push!(id, c.id) push!(name, match.name) push!(ra, match.ra) blob - 5481c545ad29ec7a14f867eaf97b348f077de8e7 blob + e6ad73652a8af14425d92f2d74bc661158576873 --- src/detection.jl +++ src/detection.jl @@ -13,6 +13,20 @@ boxes of `box_size` pixels, are extracted with `Photom detection's flux is then measured with circular aperture photometry of radius `aperture_radius` pixels on the background-subtracted image. +Each aperture is centered on a sub-pixel-refined position (a flux-weighted +centroid within `aperture_radius` of `PeakMesh`'s own integer-pixel peak +โ€” falls back to that raw peak whenever the refinement isn't trustworthy; +see [`_refine_centroid`](@ref)), not the raw integer peak itself: a real +star's true position doesn't move frame to frame, but which *integer* +pixel reads highest does, under noise โ€” and since `aperture_radius` is +comparable to a typical PSF scale, that alone changes how much flux a +fixed aperture encloses, frame to frame, for no real reason. This was the +actual cause of a real, measured false-positive floor in +[`find_variable_sources`](@ref)'s variability test (see the +[Investigation Log](https://richard7987.github.io/AsteroidPipeline.jl/dev/investigation-log#The-centroid-fix-barely-moved-the-false-positive-floor-โ€”-the-real-cause-was-a-systematic-error-floor)) โ€” the refinement only affects where the aperture +is centered, never the returned `x`/`y` (still `PeakMesh`'s own integer +position, unchanged). + Returns a table with columns `x`, `y` (pixel position), `peak` (background- subtracted peak pixel value), `flux` (aperture sum), and `flux_err`. By default (`gain=nothing`) `flux_err` is `Photometry.photometry`'s propagated @@ -57,8 +71,12 @@ function detect_sources(image::AbstractMatrix{<:Real}; # 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] + # `(row.y, row.x)` here is the fix (dim1, dim2 order); `_refine_centroid` + # takes and returns positions in that same (dim1, dim2) order. + apertures = [let (d1, d2) = _refine_centroid(subtracted, row.y, row.x, aperture_radius) + CircularAperture(d1, d2, aperture_radius) + end + 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) @@ -67,3 +85,36 @@ function detect_sources(image::AbstractMatrix{<:Real}; return Table(x=peaks.x, y=peaks.y, peak=peaks.value, flux=flux, flux_err=flux_err) end + +""" + _refine_centroid(subtracted, d1, d2, radius) -> (Float64, Float64) + +Refine an integer pixel position `(d1, d2)` (indices into `subtracted`'s +1st and 2nd dimensions respectively) to a sub-pixel flux-weighted +centroid within a `radius`-pixel window around it โ€” the same +first-moment technique [`estimate_psf`](@ref) already uses per star +stamp, applied here per detection instead. + +Falls back to the unmodified `(d1, d2)` (as `Float64`s) whenever the +refinement can't be trusted: the window would run off the array edge, +the window's total flux isn't positive, or the computed shift exceeds +`radius` itself (a sign the "centroid" is being pulled toward a neighbor +or noise excursion, not the true peak) โ€” never makes the position worse +than the raw integer peak it started from. +""" +function _refine_centroid(subtracted::AbstractMatrix{<:Real}, d1::Integer, d2::Integer, radius::Real) + r = max(1, round(Int, radius)) + n1, n2 = size(subtracted) + (r < d1 <= n1 - r && r < d2 <= n2 - r) || return Float64(d1), Float64(d2) + + stamp = subtracted[d1-r:d1+r, d2-r:d2+r] + total = sum(stamp) + total > 0 || return Float64(d1), Float64(d2) + + offs = Float64.(-r:r) + delta1 = sum(offs .* sum(stamp, dims=2)[:]) / total + delta2 = sum(offs .* sum(stamp, dims=1)[:]) / total + (abs(delta1) <= radius && abs(delta2) <= radius) || return Float64(d1), Float64(d2) + + return d1 + delta1, d2 + delta2 +end blob - 320b8e8626a390617bfaec3e4baafa730d793a7a blob + 64d8597b5d22c7e944e34c1310d9ea4f0c124c9f --- src/pipeline.jl +++ src/pipeline.jl @@ -3,7 +3,7 @@ timestamp_key::AbstractString="MJD-OBS", threshold::Real=5.0, box_size::NTuple{2,<:Integer}=(5, 5), aperture_radius::Real=3.0, max_speed::Real=Inf, - match_radius::Real=2.0, min_frames::Integer=length(fits_paths), + match_radius::Real=2.0, min_frames::Union{Nothing,Integer}=nothing, 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) @@ -66,6 +66,16 @@ 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. +A gated frame contributes zero detections, and `link_candidates` +(via `min_frames`) requires every frame to match by default โ€” so a +single real gated frame used to make no tracklet reachable at all unless +the caller manually lowered `min_frames`. `min_frames`'s default +(`nothing`) now accounts for this automatically: the value actually used +is `length(fits_paths)` minus however many frames got gated this run, +computed *after* detection, not the eager keyword default โ€” pass +`min_frames` explicitly to override this and get the old, literal +behavior. + 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 @@ -74,8 +84,8 @@ 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 +combinatorial side effect on `link_candidates`'s `min_frames` (see the +[Investigation Log](https://richard7987.github.io/AsteroidPipeline.jl/dev/investigation-log#The-quality-gate's-combinatorial-side-effect-on-tracklet-count)). 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 @@ -105,16 +115,18 @@ function run_pipeline(fits_paths::AbstractVector{<:Abs timestamp_key::AbstractString="MJD-OBS", threshold::Real=5.0, box_size::NTuple{2,<:Integer}=(5, 5), aperture_radius::Real=3.0, max_speed::Real=Inf, - match_radius::Real=2.0, min_frames::Integer=length(fits_paths), + match_radius::Real=2.0, min_frames::Union{Nothing,Integer}=nothing, 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) - detections_per_frame, wcs_per_frame, timestamps = _detect_all_frames( + detections_per_frame, wcs_per_frame, timestamps, n_gated = _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, photometric_outlier_threshold) - tracklets = link_candidates(detections_per_frame, timestamps; max_speed, match_radius, min_frames) + effective_min_frames = min_frames === nothing ? length(fits_paths) - n_gated : min_frames + tracklets = link_candidates(detections_per_frame, timestamps; + max_speed, match_radius, min_frames=effective_min_frames) return astrometric_calibrate(tracklets, wcs_per_frame, timestamps) end @@ -122,7 +134,7 @@ end _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) - -> (detections_per_frame, wcs_per_frame, timestamps) + -> (detections_per_frame, wcs_per_frame, timestamps, n_gated) The per-frame detection stage shared by [`run_pipeline`](@ref) (which links these into movers) and [`search_field`](@ref) (which additionally @@ -138,6 +150,11 @@ applied on the `reference` (ZOGY) path, since `S_corr` 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. + +`n_gated` is how many frames the `quality_max_std` gate excluded (always +`0` on the raw path, which has no such gate) โ€” both callers use it to +auto-adjust their own `min_frames` default, since a gated frame otherwise +silently makes no tracklet/variable-candidate reachable at all. """ function _detect_all_frames(fits_paths::AbstractVector{<:AbstractString}; timestamp_key::AbstractString="MJD-OBS", @@ -149,6 +166,7 @@ function _detect_all_frames(fits_paths::AbstractVector detections_per_frame = [] wcs_per_frame = WCSTransform[] timestamps = Float64[] + n_gated = 0 # Computed once, not per frame: reference.image is the same every # iteration, and this is what zogy_subtract's astrometric-noise term @@ -209,6 +227,7 @@ function _detect_all_frames(fits_paths::AbstractVector if frame_std > quality_max_std @warn "skipping frame: S_corr std exceeds quality_max_std" path frame_std quality_max_std image = permutedims(zeros(size(s_corr))) + n_gated += 1 else image = permutedims(s_corr .* valid) end @@ -237,14 +256,15 @@ function _detect_all_frames(fits_paths::AbstractVector end end - return detections_per_frame, wcs_per_frame, timestamps + return detections_per_frame, wcs_per_frame, timestamps, n_gated end """ search_field(fits_paths; , variability_position_tolerance::Real=2.0, - variability_min_frames::Integer=length(fits_paths), - variability_chi2_threshold::Real=50.0) + variability_min_frames::Union{Nothing,Integer}=nothing, + variability_chi2_threshold::Real=10.0, + variability_systematic_error_fraction::Real=0.02) -> (movers=, variables=
) Run [`run_pipeline`](@ref)'s asteroid-candidate search and @@ -257,15 +277,19 @@ since `run_pipeline` alone would need a second full pa frames to also find variables. All keywords through `photometric_outlier_threshold` are exactly -`run_pipeline`'s (see its docstring); `variability_position_tolerance`, +`run_pipeline`'s (see its docstring, including how `min_frames`'s default +now accounts for quality-gated frames); `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). +`variability_normalize`, `variability_max_relative_error`, and +`variability_systematic_error_fraction` are forwarded to +[`find_variable_sources`](@ref) as its `position_tolerance`, `min_frames`, +`chi2_threshold`, `normalize`, `max_relative_error`, and +`systematic_error_fraction`. `variability_min_frames` defaults the same +way `min_frames` does โ€” `nothing` means "every non-gated frame must +match". `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`, @@ -278,30 +302,36 @@ function search_field(fits_paths::AbstractVector{<:Abs timestamp_key::AbstractString="MJD-OBS", threshold::Real=5.0, box_size::NTuple{2,<:Integer}=(5, 5), aperture_radius::Real=3.0, max_speed::Real=Inf, - match_radius::Real=2.0, min_frames::Integer=length(fits_paths), + match_radius::Real=2.0, min_frames::Union{Nothing,Integer}=nothing, 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=50.0, + variability_min_frames::Union{Nothing,Integer}=nothing, + variability_chi2_threshold::Real=10.0, variability_normalize::Bool=(reference === nothing), - variability_max_relative_error::Real=0.10) - detections_per_frame, wcs_per_frame, timestamps = _detect_all_frames( + variability_max_relative_error::Real=0.10, + variability_systematic_error_fraction::Real=0.02) + detections_per_frame, wcs_per_frame, timestamps, n_gated = _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, photometric_outlier_threshold) - tracklets = link_candidates(detections_per_frame, timestamps; max_speed, match_radius, min_frames) + effective_min_frames = min_frames === nothing ? length(fits_paths) - n_gated : min_frames + tracklets = link_candidates(detections_per_frame, timestamps; + max_speed, match_radius, min_frames=effective_min_frames) movers = astrometric_calibrate(tracklets, wcs_per_frame, timestamps) + effective_variability_min_frames = variability_min_frames === nothing ? + length(fits_paths) - n_gated : variability_min_frames variable_groups = find_variable_sources( detections_per_frame, timestamps; position_tolerance=variability_position_tolerance, - min_frames=variability_min_frames, + min_frames=effective_variability_min_frames, chi2_threshold=variability_chi2_threshold, normalize=variability_normalize, - max_relative_error=variability_max_relative_error) + max_relative_error=variability_max_relative_error, + systematic_error_fraction=variability_systematic_error_fraction) variables = astrometric_calibrate(variable_groups, wcs_per_frame, timestamps) return (movers=movers, variables=variables) blob - 9a4a64f511e171647c73a85e950ab1dd8d7bde85 blob + 2b2b27948f0ca5c2c63b97a597cdafe8461ba09c --- src/variables.jl +++ src/variables.jl @@ -1,20 +1,48 @@ """ - variability_chi2(flux, flux_err) -> (chi2, dof) + variability_chi2(flux, flux_err; systematic_error_fraction::Real=0.02) -> (chi2, dof) -Chi-squared goodness of fit of `flux` (with per-point uncertainty -`flux_err`) against the constant-flux (non-variable) null hypothesis, -using the error-weighted mean as the constant-flux estimate. `dof` is -`length(flux) - 1` (one free parameter, the mean). +Chi-squared goodness of fit of `flux` against the constant-flux +(non-variable) null hypothesis, using the error-weighted mean as the +constant-flux estimate. `dof` is `length(flux) - 1` (one free parameter, +the mean). +The per-point error used is not `flux_err` alone, but +`sqrt(flux_err^2 + (systematic_error_fraction * flux)^2)` โ€” a systematic +error floor, proportional to flux, added in quadrature. Without it, a +bright star's tiny *formal* (Poisson/background) error makes this test +wildly oversensitive to real but non-astrophysical systematics (imperfect +flat-fielding, frame-to-frame PSF variation, ...) that never show up in +`flux_err` at all. This was found, not assumed: on real ZTF data (field +451), the false positives this chi2 test produced were *not* the +faintest, noisiest stars, as a naive read of the S/N floor would suggest +โ€” they were the **brightest** ones, with the smallest formal errors (a +median 0.25% relative error among flagged stars vs. 2.93% among the +rest), the textbook signature of a systematic floor rather than +underestimated statistical noise. A floor of 1% (a standard value in +forced-photometry pipelines, tried first) cut the false-positive rate at +a reduced-chi2 threshold of 3 from 22% (no floor) to 6%, and at +threshold 20 from 13% to essentially 0%. Sweeping the floor further +against the same real, matched stationary stars โ€” not guessed โ€” found +more room: 2% (the current default) cuts the threshold-10 rate to 2.0% +(vs 1%'s 3.3%), while a real, independently-confirmed variable +(ASASSN-V J183620.31, checked at the same floor sweep) still clears +threshold 10 with a healthy margin (reduced chi2 โ‰ˆ 31 at 2%, only +dropping below 10 once the floor reaches 5%) โ€” see +[`find_variable_sources`](@ref)'s docstring for the full before/after +table. + A large `chi2 / dof` is evidence against the null hypothesis โ€” i.e. -evidence of genuine flux variability rather than measurement noise. -Used by [`find_variable_sources`](@ref); exposed separately since a -caller may want the raw statistic rather than only a threshold decision. +evidence of genuine flux variability beyond both statistical noise and +this systematic floor. Used by [`find_variable_sources`](@ref); exposed +separately since a caller may want the raw statistic rather than only a +threshold decision. """ -function variability_chi2(flux::AbstractVector{<:Real}, flux_err::AbstractVector{<:Real}) +function variability_chi2(flux::AbstractVector{<:Real}, flux_err::AbstractVector{<:Real}; + systematic_error_fraction::Real=0.02) length(flux) == length(flux_err) || throw(ArgumentError("flux and flux_err must have the same length")) - weights = 1.0 ./ flux_err .^ 2 + eff_err = sqrt.(flux_err .^ 2 .+ (systematic_error_fraction .* flux) .^ 2) + weights = 1.0 ./ eff_err .^ 2 weighted_mean = sum(flux .* weights) / sum(weights) chi2 = sum(((flux .- weighted_mean) .^ 2) .* weights) return chi2, length(flux) - 1 @@ -99,9 +127,10 @@ end find_variable_sources(detections_per_frame, timestamps; position_tolerance::Real=2.0, min_frames::Integer=length(detections_per_frame), - chi2_threshold::Real=50.0, + chi2_threshold::Real=10.0, normalize::Bool=true, - max_relative_error::Real=0.10) + max_relative_error::Real=0.10, + systematic_error_fraction::Real=0.02) Match source detections across frames by consistent *position* (zero assumed motion) rather than [`link_candidates`](@ref)'s linear-motion @@ -115,7 +144,7 @@ matching [`link_candidates`](@ref)'s convention (used frame order meaningful for periodicity follow-up, not for any position prediction). -Two corrections are applied before any variability test: +Three corrections are applied before any variability test: - **Low-S/N floor** (`max_relative_error`, default 10%): detections with `flux_err / flux` at or above this are dropped before matching, in every @@ -128,6 +157,9 @@ Two corrections are applied before any variability tes 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). +- **Systematic error floor** (`systematic_error_fraction`, default 2%, + forwarded to [`variability_chi2`](@ref)): the actual dominant fix for + this function's real, measured false-positive floor โ€” see below. 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 @@ -136,24 +168,42 @@ frame (after its own S/N floor) is matched. Groups rea 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. +`chi2_threshold`'s default and `systematic_error_fraction` (forwarded to +[`variability_chi2`](@ref)) were both calibrated against the same real +ZTF field (451, 119 matched, high-S/N stationary stars) this whole +docstring measures against โ€” and the calibration story is worth reading, +because the first hypothesis tried here was wrong. `detect_sources` used +to center its aperture on `PeakMesh`'s raw *integer*-pixel peak; the +working theory was that per-frame pixel-grid jitter in that peak, against +a small `aperture_radius`, was producing spurious flux swings. Fixing +that (`detect_sources` now refines to a sub-pixel centroid โ€” see its own +docstring) barely moved the false-positive rate at all (23% โ†’ 22% at a +threshold of 3; 13% โ†’ 13% at 20) โ€” a real, measured non-result, not +swept under the rug. Checking *which* stars were actually being flagged +settled it: not the faintest ones, as low-S/N selection would predict, +but the **brightest** โ€” median 0.25% relative error among flagged stars +vs. 2.93% among the rest, the signature of a systematic error floor +(imperfect flat-fielding, frame-to-frame PSF variation โ€” nothing +`flux_err`'s pure Poisson/background model captures) rather than +underestimated statistical noise. Adding that floor +(`systematic_error_fraction`, in `variability_chi2`) is what actually +fixed it: a 1% floor cut it to 6% at a threshold of 3, ~0% at 20; a real, +large improvement over the 23-8% range measured before either fix, but +sweeping the floor further (against the same 119 real stars) found more +room without giving up real sensitivity โ€” checked against an actual +confirmed variable, not just the stationary-star side of the tradeoff. +`systematic_error_fraction=0.02` (the current default) cuts the +threshold-10 false-positive rate to 2.0% (vs 1%'s 3.3%), while a real, +independently-confirmed variable (ASASSN-V J183620.31 โ€” see the +[Investigation Log](https://richard7987.github.io/AsteroidPipeline.jl/dev/investigation-log#Validating-search_field-against-a-real,-independently-confirmed-variable-star)) +still clears `chi2_threshold=10.0` with a healthy margin (reduced chi2 โ‰ˆ +31, over 3x the threshold) at this floor โ€” a floor of 5% would erase that +same real signal (reduced chi2 drops to โ‰ˆ5, below threshold), so 2% is +not "as high as possible", it's the largest floor checked that still +leaves real variability comfortably detectable. Still not zero: 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. 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 @@ -170,7 +220,7 @@ What the chi2 test is discriminating against depends o 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 + mode at bright stars (see the [Investigation Log](https://richard7987.github.io/AsteroidPipeline.jl/dev/investigation-log#PSF-timing-and-astrometric-noise-bugs-behind-ZOGY's-excess-detections)), not the primary 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. @@ -192,9 +242,10 @@ 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=50.0, + chi2_threshold::Real=10.0, normalize::Bool=true, - max_relative_error::Real=0.10) + max_relative_error::Real=0.10, + systematic_error_fraction::Real=0.02) nframes = length(detections_per_frame) nframes == length(timestamps) || throw(ArgumentError("detections_per_frame and timestamps must have the same length")) @@ -221,7 +272,8 @@ function find_variable_sources(detections_per_frame, t length(group) < min_frames && continue - chi2, dof = variability_chi2([p.flux for p in group], [p.flux_err for p in group]) + chi2, dof = variability_chi2([p.flux for p in group], [p.flux_err for p in group]; + systematic_error_fraction) dof > 0 && chi2 / dof > chi2_threshold && push!(groups, group) end blob - 59a8d96fd87c5a38fe5881f1ba33b9c26c3a87f5 blob + eceb44a816f34f836fdba4e35ed55e737ccbd384 --- test/runtests.jl +++ test/runtests.jl @@ -189,6 +189,22 @@ using Reproject @test ra2 โ‰ˆ ra atol=1e-9 @test dec2 โ‰ˆ dec atol=1e-9 + # Regression test for a real bug found on real Pan-STARRS1 (IASC + # practice campaign) headers: CNPIX1/CNPIX2 (a legacy IRAF/DSS + # plate-astrometry keyword pair, unrelated to this header's real + # CTYPE/CRVAL/CRPIX/CDELT WCS) makes wcslib build a separate, + # degenerate implicit WCS and raise "Linear transformation matrix + # is singular" for the *whole* header, even though the real WCS + # parses fine alone โ€” confirmed by bisecting a real PS1 header + # down to this exact keyword pair; see load_wcs's docstring. + cnpix_cards = rpad("CNPIX1 = 0", 80) * rpad("CNPIX2 = 0", 80) + header_with_cnpix = WCS.to_header(wcs) * cnpix_cards + @test_throws "singular" WCS.from_header(header_with_cnpix) + loaded_cnpix = load_wcs(header_with_cnpix) + ra3, dec3 = pix_to_sky(loaded_cnpix, 500.0, 500.0) + @test ra3 โ‰ˆ ra atol=1e-9 + @test dec3 โ‰ˆ dec atol=1e-9 + tracklets = [[(frame=1, x=500.0, y=500.0), (frame=2, x=501.0, y=500.0)]] timestamps = [2460000.5, 2460000.51] candidates = astrometric_calibrate(tracklets, [wcs, wcs], timestamps) @@ -557,6 +573,18 @@ using Reproject @test !(3 in candidates.frame) # the glow-patch frame contributed nothing @test 1 in candidates.frame # the real transient still recovered elsewhere + # same real transient, but with min_frames *omitted* entirely + # (its new default, nothing, auto-subtracts the 1 gated frame + # from length(scipaths)=3, giving an effective min_frames=2) โ€” + # must reach the same real result as the explicit min_frames=1 + # above, without the caller having to know a frame was gated. + auto_candidates = @test_logs (:warn, r"quality_max_std") match_mode = :any run_pipeline( + scipaths; threshold=6.0, match_radius=5.0, + reference=reference, psf_threshold=20.0, psf_min_separation=15.0) + @test !isempty(auto_candidates) + @test !(3 in auto_candidates.frame) + @test 1 in auto_candidates.frame + # lowering the threshold further should also drop the clean frames @test_logs((:warn, r"quality_max_std"), (:warn, r"quality_max_std"), (:warn, r"quality_max_std"), match_mode = :any, @@ -775,6 +803,27 @@ using Reproject end try + # regression test for batched (OR-chained, one request for all + # candidates) cross-matching: 3 real, independently-verified + # positions in one call โ€” 2 real VSX variables (same RS star + # as above, plus V0651 Ori, type EW) and 1 real non-match (near + # the celestial pole) โ€” confirms each returned row is + # attributed to the *correct* candidate id from a single + # shared request, not just that matches exist somewhere. + batch = [(id=10, ra=36.2344, dec=2.06997), + (id=20, ra=83.1937, dec=5.41603), + (id=30, ra=0.0, dec=89.9)] + matches = crossmatch_catalog(batch, :vsx; radius=5.0) + @test length(matches) == 2 + @test Set(matches.id) == Set([10, 20]) + @test only(matches[matches.id.==10].class) == "RS" + @test only(matches[matches.id.==20].class) == "EW" + 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.