Guides

Fill known bad pixels

gspice.infill() replaces known bad pixels or user-selected wavelength intervals with their conditional Gaussian prediction from the remaining good pixels of each spectrum.

Set up

The examples share one model, fitted to 2000 simulated spectra of 100 pixels. They predict pixels of 500 other spectra. The Quickstart explains each call.

import numpy as np

import gspice

rng = np.random.default_rng(146)
train = gspice.simulate(2000, 100, rng=rng)
test = gspice.simulate(500, 100, rng=rng)
model, _ = gspice.fit(train.flux, train.ivar)

The 500 test spectra are held-out spectra: the model was not fitted to them.

Replace the bad pixels

Suppose that your pipeline marks some pixels of each spectrum as bad. A bad pixel has a nonzero mask or ivar == 0. Here a random 2% of the pixels are bad:

bad = rng.random(test.flux.shape) < 0.02
result = gspice.infill(model, test.flux, test.ivar, mask=bad)
print(f"{result.replaced.sum()} pixels replaced")
print(np.array_equal(result.flux[~bad], test.flux[~bad]))

1033 pixels replaced
True

infill replaced every bad pixel by its prediction from the good pixels of its spectrum. The other pixels keep their input values. What you get back lists every field of result.

Predict a range of wavelengths

Sometimes you need a prediction where the pixels are good, for example to compare it with an absorption line. Pass the pixels to predict as replace:

gap = np.zeros(100, dtype=bool)
gap[52:72] = True
result = gspice.infill(model, test.flux, test.ivar, replace=gap)
print(result.replaced.sum(axis=1)[:5])

[20 20 20 20 20]

A replace mask of shape (npix,) applies to every spectrum. Here it holds pixels 52 to 71, around the absorption line at pixel 62. With a wavelength array wave, the mask (wave > 4840) & (wave < 4890) selects a range in Å.

The prediction uses only the conditioning pixels, the pixels outside replace that are good. A bad pixel outside replace keeps its input value, and the prediction does not use it.

Get the uncertainty

The variance of each pixel

A prediction is uncertain. result.var holds the variance of each replaced pixel, in flux², and includes the noise of the pixel. The simulated pixels in the gap are good, so you can compare their observed flux with the prediction:

pull = (test.flux - result.flux)[:, gap] / np.sqrt(result.var[:, gap])
print(f"standard deviation of the pulls: {np.std(pull):.2f}")

standard deviation of the pulls: 1.03

A pull is the observed flux minus the prediction, divided by the square root of the variance. A standard deviation above 1 means that the variances are too small.

The option inflate=True multiplies the variances by the finite-sample factor, by which the errors of the model widen the residuals of held-out spectra:

result = gspice.infill(model, test.flux, test.ivar, replace=gap, inflate=True)
pull = (test.flux - result.flux)[:, gap] / np.sqrt(result.var[:, gap])
print(f"standard deviation of the pulls: {np.std(pull):.2f}")

standard deviation of the pulls: 0.98

The finite-sample factor explains the factor, and which models support inflate.

The covariance of the pixels

The errors of neighboring predicted pixels are correlated. A measurement that combines several pixels, such as the total flux in a band, needs these correlations. With covariance=True, result.covariance(s) returns the replaced pixels of spectrum s and their covariance $C$ :

result = gspice.infill(
    model, test.flux, test.ivar, replace=gap, inflate=True, covariance=True
)
pixels, cov = result.covariance(0)
w = np.ones(pixels.size)
band = w @ result.flux[0, pixels]
print(f"band flux {band:.3f} ± {np.sqrt(w @ cov @ w):.3f}")
print(f"without the correlations: ± {np.sqrt(w**2 @ np.diag(cov)):.3f}")

band flux 18.966 ± 0.243
without the correlations: ± 0.133

The band flux is $w^\mathsf{T} f$ , with the weights $w$ and the predicted flux $f$ . Its variance is $w^\mathsf{T} C w$ . The diagonal of $C$ alone gives a smaller uncertainty.

Over all 500 test spectra, the covariance gives the right spread:

z, z_diagonal = [], []
for s in range(len(test.flux)):
    pixels, cov = result.covariance(s)
    error = np.sum(test.flux[s, pixels] - result.flux[s, pixels])
    z.append(error / np.sqrt(cov.sum()))
    z_diagonal.append(error / np.sqrt(np.trace(cov)))
print(f"{np.std(z):.2f} with the covariance, {np.std(z_diagonal):.2f} without")

1.01 with the covariance, 1.84 without

Each value is the standard deviation of the band-flux error over the 500 spectra, in units of its predicted uncertainty.

Use draws

result.flux holds the conditional mean: the average of all spectra that agree with the conditioning pixels. An average is smoother than any one spectrum. A draw is one plausible spectrum instead. Pass a numpy.random.Generator as rng to get one in result.draws:

result = gspice.infill(model, test.flux, test.ivar, replace=gap, rng=rng)
observed = np.std(test.flux[:, gap], axis=0)
for name, flux in [("mean", result.flux), ("one draw", result.draws)]:
    spread = np.std(flux[:, gap], axis=0) / observed
    print(f"{name:>8}: {spread.mean():.2f}")

    mean: 0.50
one draw: 1.03

Each value compares the spread of the flux over the 500 test spectra with that of the observed flux, at each pixel of the gap, and averages the ratio over the gap. The mean has about half the spread of real spectra. A draw has about the same spread.

So choose by the next step:

Plot a repaired spectrum, or make a measurement linear in the flux: use the mean. covariance=True gives the uncertainty of the measurement exactly, $w^\mathsf{T} C w$ .
Treat the spectrum as an observation: use one draw. Examples are a χ² fit that assumes the noise, a noise estimate, or a classifier trained on real spectra. The smooth mean would bias these steps.
Get the uncertainty of a nonlinear measurement: use several draws. Run infill on copies of the spectrum, measure each draw, and take the spread.

The last method is called multiple imputation. This example measures the depth of the line in the gap of spectrum 0, from 20 draws:

copies = np.repeat(test.flux[:1], 20, axis=0)
copies_ivar = np.repeat(test.ivar[:1], 20, axis=0)
result = gspice.infill(model, copies, copies_ivar, replace=gap, rng=rng)
draws = list(result.draws)
depths = [1 - draw[gap].min() for draw in draws]
print(f"line depth {np.mean(depths):.3f} ± {np.std(depths):.3f}")

line depth 0.244 ± 0.040

The depth is 1 minus the lowest flux in the gap. The spread of the 20 depths is the uncertainty that the missing pixels add.

0.40.60.81.00255075100

observed flux
masked
prediction
draw
flux without noise

The conditional mean and 3 draws of the multiple imputation replace pixels 52 to 71 of test.flux[0], at a signal-to-noise ratio of 44. The mean is smooth. Each draw scatters from pixel to pixel like the observed flux outside the gap.