commit - dadb629afdb72799c3af862cb777266bede2d495
commit + 8ff5933c85c68fcb18d5be2ee05c05d4d2c6df00
blob - 014e8d58f76ef8cc3b461a0a3d1a9c34794a23f2
blob + 413f5933626675d76786f401e56e3c9361681dbf
--- docs/src/index.md
+++ docs/src/index.md
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
`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
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
`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
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
`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
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
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
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
"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
+#=
+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/<set name>/*.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
+#=
+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
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
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)
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
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)
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
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[]
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
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
# 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)
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
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)
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
`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
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
_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
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",
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
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
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; <all run_pipeline keywords>,
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=<table>, variables=<table>)
Run [`run_pipeline`](@ref)'s asteroid-candidate search and
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`,
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
"""
- 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
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
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
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
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
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.
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"))
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
@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)
@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,
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.
