fixedbugs

Find and mask the fixed bugs in APERO t.fits spectra.

Documentation with figures: https://eartigau.github.io/fixedbugs/

---

What this is

APERO, the reduction pipeline of SPIRou and NIRPS, writes one *t.fits file per exposure. Inside it, the science flux extension holds the extracted, blaze-shaped, telluric corrected 2D spectrum: one row per diffraction order, one column per detector pixel.

That spectrum is never perfect. A handful of pixels in every order carry a residual that no amount of averaging removes, because the residual sits at the same detector pixel in every exposure. The usual culprits:

• imperfectly corrected telluric absorption (deep water or methane bands),
• residuals of the OH airglow emission subtraction,
• persistence left over from the daytime calibrations, above all by the saturated lines of the Uranium-Neon wavelength-calibration lamp, which stay imprinted on the detector well into the night,
• other detector features: bad pixel clusters, ghosts, order edges.

All of these are fixed in the observer's rest frame. The star is not: over a year the barycentric Earth radial velocity (BERV) sweeps the stellar spectrum back and forth by up to ±30 km/s, which at SPIRou's 2.28 km/s per pixel is a couple of dozen pixels. Stack enough exposures and the two separate cleanly. That is the whole trick, and it is why we call these residuals fixed bugs: they are bugs, and they never move.

fixedbugs measures where they are and sets those pixels to NaN, writing a clean copy of every input file. Downstream tools that already skip NaNs (LBL, CCF codes, template builders) then ignore them for free.

321 SPIRou exposures of TOI-2120, the middle 1000 pixels of one order. Each row is one exposure. The vertical stripes are pixels that misbehave in every single exposure: they do not move with the star, so they are not the star. The bar underneath is the mask fixedbugs derives from them.

How it works

Two passes over the data.

Pass 1, build the mask. For every exposure and every order:

1. Shift a high signal-to-noise 1D template of the same star (an APERO s1d_v template, which lives in the barycentric rest frame) into the observer frame using that exposure's BERV, and multiply it by the blaze. That is the model of what the exposure should look like.
2. Divide the observed flux by that model, normalise by the median. A perfect reduction would give a flat line at 1.
3. Subtract a running median, to throw away the large-scale mismatch we do not care about: template colour, imperfect blaze, continuum normalisation.
4. Divide by a robust sigma from the 16th–84th percentile spread, so every exposure contributes on the same scale whatever its signal-to-noise.

Then take, pixel by pixel, the median of the absolute value across all exposures. A pixel that misbehaves once (a cosmic ray) is voted down by the median; a pixel that misbehaves every time survives.

That curve is not flat even when nothing is wrong, because the noise itself grows towards the order edges where the blaze drops. So rather than one global threshold, fit_lower_envelope fits the noise floor along the order: a degree-2 polynomial under an asymmetric loss, where a point above the curve costs asymmetry times its squared residual and a point below costs 1 - asymmetry times it. At the default asymmetry: 0.02 that is a 49:1 ratio, so the curve settles just under the bulk of the points instead of running through their middle, and the bad pixels can no longer drag their own reference level upwards and hide behind it.

The small penalty on positive residuals is deliberately not zero. If it were, nothing would stop the fit sliding down forever, since going lower would cost nothing. Keeping a little cost per point means a curve that has sunk below everything pays for every point at once, which pulls it back up to the floor. That single term is what makes the problem well posed. The fit is solved by iteratively reweighted least squares and converges in 5 to 10 passes; it is implemented in fixedbugs.py itself, so there is no fitting library to install.

Pixels exceeding that floor by more than nsig_mask are flagged, so the default nsig_mask: 2.0 means "more than three times the local noise floor". Flagged regions are then grown by dilation pixels, because a bad pixel bleeds into its neighbours through the extraction profile.

One order. Grey is the median absolute residual over all 321 exposures; red is the fitted noise floor, hugging the bottom of the point cloud and following its natural rise towards the red end; pink bands are the pixels flagged as being in excess of it.

Pass 2, apply the mask. The mask is stored as an array of 0.0 (keep) and NaN (reject) and simply added to the flux: adding zero changes nothing, adding NaN gives NaN. Every other extension and the full header are copied verbatim, so the output is still a valid APERO product.

SPIRou and NIRPS

The two instruments do not name their science fibre the same way:

Instrument	Science fibre	Extensions	Typical shape
SPIRou	AB (polarimetric fibres A and B combined)	FluxAB, WaveAB, BlazeAB	49 × 4088
NIRPS	A (B looks at the sky)	FluxA, WaveA, BlazeA	75 × 4088

fixedbugs reads the INSTRUME keyword of the first input file and picks the right one on its own. Nothing to change when you switch instrument. You can still force a fibre (extensions.fiber: "C", say) if you want to clean a calibration or sky channel.

Array dimensions are always measured from the data, never assumed, and the speed of light comes from astropy.constants rather than a retyped literal (BERV is in km/s in APERO headers).

Install

fixedbugs is a single script. Clone it, install its dependencies, and run it:

git clone https://github.com/eartigau/fixedbugs.git
cd fixedbugs
python -m venv venv && source venv/bin/activate   # keeps this off your system Python
pip install -r requirements.txt

requirements.txt pulls in numpy, scipy, astropy, matplotlib, pyyaml and tqdm, and nothing else: the polynomial fit that sets the detection threshold lives in fixedbugs.py.

Use

Everything lives in a YAML file. Copy config.yaml, point the three required entries at your own data, and run:

cp config.yaml my_star.yaml
# edit my_star.yaml: io.input_glob, io.output_dir, template.file
python fixedbugs.py my_star.yaml

Paths in the YAML may be absolute, or relative to the YAML file itself (not to wherever you launch the script from), so a configuration keeps working from any directory.

You need three things:

• io.input_glob — a shell pattern matching the t.fits files. They must all be the same star on the same instrument, so that one template and one detector geometry apply.
• io.output_dir — where the masked copies go. One output per input, same filename.
• template.file — the APERO s1d_v template of that star (Template_s1dv_<OBJECT>_sc1d_v_file_AB.fits on SPIRou, ..._file_A.fits on NIRPS).

Useful flags:

Flag	What it does
--mask-only	build the mask and the figures, write no masked files
--apply-only	reuse an existing mask instead of measuring a new one
--mask FILE	which mask to reuse with --apply-only
--no-plots	skip the diagnostic figures

--apply-only is how you carry a mask built on a long, well sampled time series over to a handful of new exposures of the same star, where there would not be enough BERV coverage to measure one.

Progress is reported as timestamped, colour-coded lines: green for what is happening, blue for a measured number, orange for something skipped, red for why a run stopped.

260904 16:58:00.16 | 321 file(s) matching /.../data/tfits/*t.fits
260904 16:58:00.16 | instrument SPIRou detected, using science fibre AB
260904 16:58:00.24 | data shape measured from 2811170t.fits: 49 orders x 4088 pixels
260904 16:58:00.24 | residual cube: 0.26 GB in memory (float32)
260904 16:58:48.24 | 50846 of 200312 pixels were already unmeasurable (25.4%: NaN in the input, or outside the template coverage)
260904 16:58:48.24 | 26170 of 149466 measurable pixels masked (17.51%)
260904 16:58:48.24 | masked fraction per order: median 12.7%, min 4.8%, max 56.9%

What comes out

The masked files, one per input, in io.output_dir. Identical to the input except that the science flux extension has NaNs where the fixed bugs were. Provenance is recorded in the header: FBVERS, FBDATE, FBINSTR, FBFIBER, FBNFILE, FBNSIG, FBPORD, FBASYM, FBNPIX, FBFRAC and the rest, so a file found on disk a year from now can be traced back to the settings that produced it.

The mask, a multi-extension FITS file:

Extension	Contents
MASK	the mask itself: 0.0 keep, NaN reject
OFFENDER	median |residual| over all exposures, in sigma
EXPECTED	the fitted noise floor: the lower envelope of OFFENDER
EXCESS	OFFENDER / EXPECTED - 1, the quantity thresholded

Keeping the three diagnostics next to the mask means you can revisit the threshold without redoing pass 1.

The figures, in io.plot_dir (PDF by default, plots.format: "png" for bitmaps):

• mask_map — the whole mask, one row per order.
• order_NN — the offender curve, the fitted expectation, the flagged pixels. This is the plot to look at when tuning nsig_mask.
• river_NN — every exposure stacked for one order, zoomed on the middle 1000 pixels (plots.river_half_width, null for the whole order). The stripes.
• before_after — one real order before and after masking.

Tuning

Two knobs matter; the rest you can leave alone.

Every percentage below counts only the measurable pixels. About a quarter of a SPIRou t.fits is already NaN before fixedbugs touches it (deep telluric bands, pixels APERO rejected, wavelengths the template does not cover), and those were never ours to reject. On the 321-exposure data set: 50 846 of the 200 312 pixels (25.4%) are unmeasurable, leaving 149 466 that fixedbugs actually judges.

detection.nsig_mask sets how aggressive the masking is. It is a relative excess over the fitted noise floor, so 2.0 means "more than three times the local noise floor". At the default dilation: 2:

nsig_mask	measurable pixels masked
0.5	49.4%
1.0	30.6%
2.0 (default)	17.5%
3.0	11.6%
4.0	7.9%
6.0	3.9%

detection.dilation is how far each flagged region is grown, in pixels. It costs more than people expect, because most flagged regions are only one or two pixels wide to begin with. At nsig_mask: 2.0:

dilation	measurable pixels masked
0	8.8%
1	13.5%
2 (default)	17.5%
5	27.3%
10	40.0%

Setting dilation: 0 leaves the flagged pixels exactly as they are; note that scipy reads iterations <= 0 as "repeat until nothing changes", so fixedbugs skips the call entirely rather than passing zero through.

detection.envelope_asymmetry controls how far under the points the noise floor sits. Measured as the share of pixels ending up below the fitted curve:

envelope_asymmetry	pixels below the curve
0.005	6.5%
0.02 (default)	15.8%
0.10	36.0%
0.40	60.5%

It must stay strictly above 0 and below 0.5; 0.5 is the symmetric case, an ordinary least-squares fit.

If an order cannot be measured at all (fewer than min_valid_fraction of its pixels are finite), fixedbugs leaves it alone and says so in the log, rather than rejecting it. Two very different things land in that bucket: a saturated telluric band, where the data really is unusable, and a template that does not reach that far in wavelength, where the data may be perfectly fine and we simply cannot check it. Since they are indistinguishable from the residuals alone, the default is to only ever remove what there is evidence against. Set detection.unusable_order_action: "mask" to reject them instead.

Memory and runtime

Pass 1 holds the whole residual cube in memory: n_files × n_orders × n_pixels × 4 bytes with detection.float32: true. 321 SPIRou exposures (49 × 4088) is 0.26 GB and takes about 50 seconds on a laptop. 1000 NIRPS exposures (75 × 4088) would be 1.2 GB. If that is too much, detection.max_files measures the mask on a subset, and --apply-only then applies it to everything.

A note for students

Four ideas in this script are worth stealing for your own work.

Use the thing that moves to find the thing that does not. Nothing here knows what a telluric line or a bad pixel looks like. The separation comes entirely from the fact that one signal is fixed in the detector frame and the other is not, and that the Earth's orbit does the sweeping for free.

Take the median of absolute values, not the mean. Averaging residuals lets one bad exposure dominate. The median over exposures asks a different question: "is this pixel misbehaving most of the time?" A cosmic ray cannot answer yes.

Fit the floor, not the middle. The residual level rises towards the order edges where the blaze drops, so a single flat cut would flag whole edges while missing bad pixels in the bright middle. Worse, an ordinary fit of that level is itself pulled upwards by the very pixels we are hunting, which then hide behind their own contribution. An asymmetric loss that fits the lower envelope breaks that circularity in a handful of lines.

Quote the honest denominator. A quarter of these files is already NaN before we start, so "13.1% of the detector masked" quietly credits fixedbugs with pixels it never touched. What a reader actually wants to know is what fraction of the usable data was removed, which is 17.5%. Pick the denominator that answers the question, and say which one you picked.

License

MIT, see LICENSE.