Legacy interface

gspice.legacy preserves the 2022 functions with the arithmetic of commit 1fbab9f so existing code continues to run during migration. Migrate from the 2022 package gives their replacements in the new API. The 2022 modules also remain available as gspice.legacy.specUtils, gspice.legacy.matrixUtils and gspice.legacy.djs_maskinterp, so each function below can be imported from its 2022 module path.

Scaling and covariance

gspice.legacy.standard_scale

Python

gspice.legacy.standard_scale(spec, ivar, mask=None)

Interpolate over bad pixels, scale spectra to about unit noise, and center.

Dvec = refscale * maskinterp(spec) - refmean, where maskinterp interpolates along each spectrum, refscale is the square root of the mean ivar of the good pixels of each spectrum, and refmean is the mean of the scaled spectra.

Parameters:

• spec (ndarray, shape (nspec, npix)) – Flux. An integer dtype truncates the interpolated values.
• ivar (ndarray, shape (nspec, npix)) – Inverse variance of the flux.
• mask (ndarray, shape (nspec, npix), optional) – Nonzero, NaN included, marks a bad pixel. None marks the pixels where ivar == 0.

Returns:

• Dvec (ndarray, shape (nspec, npix)) – The scaled and centered spectra, dimensionless.
• refscale (ndarray, shape (nspec,)) – The scale of each spectrum, in units of 1/flux. A non-finite ivar at a good pixel makes the scale of its spectrum non-finite, and with it every value of Dvec and refmean.
• refmean (ndarray, shape (npix,)) – The mean of the scaled spectra.

Raises:

• ValueError – If a spectrum has no good pixel.

gspice.legacy.covar

Python

gspice.legacy.covar(spec, checkmean=False)

Sample covariance of the pixels over the spectra.

cov = (spec - refmean)^T (spec - refmean) / (nspec - 1), where refmean is the mean spectrum.

Parameters:

• spec (ndarray, shape (nspec, npix)) – Spectra, float64.
• checkmean (bool, optional) – If checkmean is True, print the minimum, maximum and standard deviation of the column means after centering.

Returns:

• cov (ndarray, shape (npix, npix)) – The covariance, in the units of spec squared.
• refmean (ndarray, shape (npix,)) – The mean spectrum, in the units of spec.

Raises:

• ValueError – If spec is not float64.

gspice.legacy.covar_iter_mask

Python

gspice.legacy.covar_iter_mask(
    flux,
    ivar,
    mask,
    nsigma=np.array([20, 8, 6]),
    maxbadpix=64,
)

Covariance of the spectra after passes of get_chimask, one per threshold.

Spectra with more than maxbadpix bad pixels in mask are dropped. Each pass runs get_chimask with the bad pixels of mask and the flags of the previous pass masked, not those of earlier passes. finalmask combines the mask of the last pass with its flags, and the covariance comes from standard_scale and covar with finalmask.

Parameters:

• flux (ndarray, shape (nspec, npix)) – Flux.
• ivar (ndarray, shape (nspec, npix)) – Inverse variance of the flux.
• mask (ndarray, shape (nspec, npix), or None) – Nonzero, NaN included, marks a bad pixel. None marks the pixels where ivar == 0.
• nsigma (sequence of float, optional) – The thresholds of the passes, in order; at least one.
• maxbadpix (int, optional) – The largest number of bad pixels in mask that a kept spectrum has.

Returns:

• cov (ndarray, shape (npix, npix)) – The covariance of the kept spectra, scaled as by standard_scale.
• finalmask (ndarray of float64, shape (nspec, npix)) – 1.0 for a masked pixel and 0.0 for another. Every pixel of a dropped spectrum is 1.0.

Raises:

• ValueError – If fewer than npix spectra are kept, or if mask, the mask of a pass or finalmask leaves a kept spectrum without a good pixel.
• numpy.linalg.LinAlgError – If the covariance of a pass is not positive definite. With exactly npix kept spectra it is singular, so a pass usually raises this error, and otherwise the returned covariance is singular.

Notes

Cost: that of get_chimask per pass, plus standard_scale and covar.

Prediction

gspice.legacy.get_chimask

Python

gspice.legacy.get_chimask(flux, ivar, mask, nsigma)

Flag the pixels that deviate from their GSPICE prediction by nsigma or more.

Z = (Dvec - pred) / sqrt(predvar), where Dvec comes from standard_scale, and pred and predvar from gp_interp at nguard 20 with the covariance of Dvec. A pixel with |Z| >= nsigma is flagged, together with its neighbors in its spectrum.

Parameters:

• flux (ndarray, shape (nspec, npix)) – Flux.
• ivar (ndarray, shape (nspec, npix)) – Inverse variance of the flux.
• mask (ndarray, shape (nspec, npix), or None) – Bad pixels, as in standard_scale.
• nsigma (float) – Threshold on |Z|.

Returns:

• chimask (ndarray of bool, shape (nspec, npix)) – True for a flagged pixel.

Raises:

• ValueError – If a spectrum has no good pixel.

Notes

Cost: that of covar and of gp_interp at nguard 20.

gspice.legacy.gp_interp

Python

gspice.legacy.gp_interp(spec, cov, nguard=20, irange=None, bruteforce=False)

Predict each pixel from the pixels outside a guard window around it (GSPICE).

For pixel i, gaussian_estimate conditions on every pixel outside i - nguard to i + nguard.

Parameters:

• spec (ndarray, shape (nspec, npix)) – Spectra, such as the output of standard_scale.
• cov (ndarray, shape (npix, npix)) – Positive definite covariance of the pixels.
• nguard (int, optional) – Half width of the guard window.
• irange (sequence of int, optional) – Predict the pixels irange[0] to irange[1], both included, and ignore later elements. The default is every pixel.
• bruteforce (bool, optional) – If True, gaussian_estimate inverts the covariance of the conditioning pixels of each predicted pixel with cholesky_inv, instead of updating the inverse of cov.

Returns:

• pred (ndarray, shape (nspec, n)) – The predictions of the n = irange[1] - irange[0] + 1 predicted pixels, in the units of spec.
• predvar (ndarray, shape (nspec, n)) – The conditional variances, in the units of spec squared. All rows are equal.

Raises:

• ValueError – If cov is not square, does not match spec or has a non-finite entry, spec has no pixels, nguard is negative, or irange is not 1-D, has fewer than 2 elements or breaks 0 <= irange[0] <= irange[1] <= npix - 1.
• TypeError – If nguard, irange[0] or irange[1] is not an integer, or is a numpy bool.
• numpy.linalg.LinAlgError – If cov is not positive definite.

Notes

Cost: about npix^3 flops for the inverse of cov, and per predicted pixel about 2 npix^2 + w^3 + 2 w (w + npix) + 4 nspec npix flops, where w is the number of pixels in the guard window, at most 2 nguard + 1. With bruteforce, per predicted pixel about (npix - w)^3 + 2 nspec npix flops.

gspice.legacy.gaussian_estimate

Python

gspice.legacy.gaussian_estimate(icond, ipred, cov, Dvec, covinv, bruteforce=False)

Gaussian conditional mean and covariance of some pixels, given others.

With S = cov, the conditioning pixels k and the predicted pixels k*, the mean is Dvec[:, k] S_kk^-1 S_kk* and the covariance is S_k*k* - S_k*k S_kk^-1 S_kk*. The default path updates covinv with submatrix_inv_mult; bruteforce inverts S_kk with cholesky_inv. If icond and ipred overlap, each shared pixel gets its own value as the mean and a variance of 0 to rounding, which can be negative.

Parameters:

• icond (ndarray, shape (npix,)) – Nonzero for the conditioning pixels k.
• ipred (ndarray, shape (npix,)) – Nonzero for the predicted pixels k*.
• cov (ndarray, shape (npix, npix)) – Covariance of the pixels.
• Dvec (ndarray, shape (nspec, npix)) – Spectra, such as the output of standard_scale.
• covinv (ndarray, shape (npix, npix)) – Inverse of cov. bruteforce does not read it.
• bruteforce (bool, optional) – If True, return (predcovar, predkstar, kstar) only.

Returns:

• predoverD (ndarray, shape (npix, nk*)) – S_kk^-1 S_kk* in the rows k and 0 in the other rows, so that the mean is Dvec @ predoverD.
• predcovar (ndarray, shape (nk*, nk*)) – The conditional covariance, in the units of cov.
• predkstar (ndarray, shape (nspec, nk*)) – The conditional mean, in the units of Dvec. With bruteforce and nk* = 1, its shape is (nspec,).
• kstar (ndarray of int, shape (nk*,)) – The indices of the predicted pixels.

Notes

Cost: about 2 npix (npix + nk*) nk* + nr^3 + 2 nr (nr + npix) nk* + 2 nspec npix nk* flops, where nr is the number of pixels outside k; with bruteforce, about nk^3 + 2 nk nk* (2 nk + nk*) + 2 nspec n nk* flops, where n is npix for one predicted pixel and nk for more.

Matrix utilities

gspice.legacy.cholesky_inv

Python

gspice.legacy.cholesky_inv(M)

Inverse of a positive definite matrix, from its Cholesky factor.

M = L L^T with L lower triangular, and LAPACK dpotri forms M^-1 from L.

Parameters:

• M (array_like, shape (N, N)) – Positive definite matrix with finite entries. Only its lower triangle enters the result.

Returns:

• ndarray of float64, shape (N, N) – The inverse of M, symmetric. A float32 M is factorized in float32.

Raises:

• ValueError – If M is not square or has a non-finite entry.
• numpy.linalg.LinAlgError – If M is not positive definite.

Notes

Cost: about N^3 flops.

gspice.legacy.submatrix_inv

Python

gspice.legacy.submatrix_inv(M, Minv, imask, bruteforce=False)

Inverse of the submatrix of M that imask keeps, from the inverse of M.

With the kept indices k and the removed indices r, A = M[k, k] and A^-1 = P - Q U^-1 Q^T, where P = Minv[k, k], Q = Minv[k, r] and U = Minv[r, r].

Parameters:

• M (ndarray, shape (N, N)) – Symmetric positive definite matrix.
• Minv (ndarray, shape (N, N)) – Inverse of M. bruteforce does not read it.
• imask (ndarray, shape (N,)) – Nonzero for the rows and columns to keep.
• bruteforce (bool, optional) – If True, invert A with cholesky_inv.

Returns:

• Ainv (ndarray, shape (nk, nk)) – The inverse of A.

Raises:

• ValueError – If M is not square or imask.shape[0] is not N.
• numpy.linalg.LinAlgError – If U is singular, or with bruteforce, if A is not positive definite.

Notes

Cost: about 2 nr^3 + 2 nk nr (nk + nr) flops, for nr removed and nk kept indices; with bruteforce, about nk^3 flops.

gspice.legacy.submatrix_inv_mult

Python

gspice.legacy.submatrix_inv_mult(M, Minv, imask, Y, MinvY, pad=True, bruteforce=False)

Product of the inverse of a submatrix of M with Y, from the inverse of M.

With the kept indices k, the removed indices r, A = M[k, k], P = Minv[k, k], Q = Minv[k, r] and U = Minv[r, r], A^-1 Y[k] = P Y[k] - Q U^-1 Q^T Y[k]. As MinvY[k] = P Y[k] + Q Y[r] and MinvY[r] = Q^T Y[k] + U Y[r], this is MinvY[k] - Q (U^-1 (MinvY[r] - U Y[r]) + Y[r]), which the function evaluates, so Y[r] can hold any values.

Parameters:

• M (ndarray, shape (N, N)) – Symmetric positive definite matrix.
• Minv (ndarray, shape (N, N)) – Inverse of M. bruteforce does not read it.
• imask (ndarray, shape (N,)) – Nonzero for the rows and columns to keep.
• Y (ndarray, shape (N, m)) – The matrix to multiply.
• MinvY (ndarray, shape (N, m)) – Minv @ Y. bruteforce checks only that it is 2-D.
• pad (bool, optional) – If True, return N rows, with 0 in the rows r.
• bruteforce (bool, optional) – If True, invert A with cholesky_inv, and return nk rows whatever pad is.

Returns:

• ndarray, shape (N, m) if pad, else (nk, m) – A^-1 Y[k].

Raises:

• ValueError – If M is not square, imask.shape[0] is not N, or Y or MinvY is not 2-D.
• numpy.linalg.LinAlgError – If U has more than one row and is not positive definite, or with bruteforce, if A is not positive definite.

Notes

Cost: about nr^3 + 2 nr (2 nr + nk) m flops, for nr removed and nk kept indices; with bruteforce, about nk^3 + 2 nk^2 m flops.

Interpolation

gspice.legacy.maskinterp1

Python

gspice.legacy.maskinterp1(yval, mask)

Interpolate linearly over the bad entries of a 1-D array.

The grid is regular. Entries before the first or after the last good entry get the value of that entry, as with the /const option of the IDL routine. With one good entry, of value v, a bad entry with a non-finite value gets v, and every other entry gets yval * 0 + v.

Parameters:

• yval (array_like, shape (n,)) – Values. An integer dtype truncates the interpolated values.
• mask (ndarray, shape (n,)) – Nonzero after conversion to int marks a bad entry, so 0.5 is good.

Returns:

• ndarray, shape (n,) – A copy of yval, interpolated. Without bad or without good entries, the copy is unchanged.

gspice.legacy.maskinterp

Python

gspice.legacy.maskinterp(yval, mask, axis)

Apply maskinterp1 to each row or to each column of a 2-D array.

Parameters:

• yval (ndarray, shape (nrow, ncol)) – Values.
• mask (ndarray, shape (nrow, ncol)) – Bad entries, as in maskinterp1.
• axis ({0, 1}) – 0 interpolates along each row, and 1 along each column, as in the IDL routine.

Returns:

• ndarray, shape (nrow, ncol) – A copy of yval, interpolated.

Raises:

• ValueError – If yval or mask is not 2-D, their shapes differ, or axis is not 0 or 1.