ExtractionEngine

class jwst.extract_1d.soss_extract.atoca.ExtractionEngine(wave_map, trace_profile, throughput, kernels, wave_grid, mask_trace_profile, global_mask=None, orders=None, threshold=0.001)

Bases: object

Run the ATOCA algorithm.

This class is basically a fitter. On instantiation, it generates a model of the detector, including a mapping between the detector pixels and the wavelength for each spectral order, the throughput and convolution kernel, and known detector bad pixels. This does not require any real data. When called, it ingests data and associated errors, then generates an output 1-D spectrum that explains the pixel brightnesses in the data within the constraints of the model.

The engine can also run in reverse: The rebuild() method generates a synthetic 2-D detector ‘observation’ from a known or fitted spectrum, and the compute_likelihood() method compares the synthetic data to the real data to generate a likelihood. This allows for a likelihood-based optimization of the spectrum.

This version models the pixels of the detector using an oversampled trapezoidal integration.

Parameters:

wave_maplist or ndarray

A list or array of 2-D arrays of the central wavelength position for each order on the detector. Has shape (N_ord, N, M). It has to have the same (N, M) as data.

trace_profilelist or ndarray

A list or array of 2-D arrays of the spatial profile for each order Has shape (N_ord, N, M). on the detector. It has to have the same (N, M) as data.

throughputlist of array or callable

A list of functions or array of the throughput at each order. If callable, the functions depend on the wavelength. If array, projected on wave_grid. Has shape (N_ord [, N_k]).

kernelscallable, sparse matrix, or None

Convolution kernel to be applied on spectrum (f_k) for each orders. Can be a callable with the form f(x, x0) where x0 is the position of the center of the kernel. In this case, it must return a 1D array (len(x)), so a kernel value for each pairs of (x, x0). If callable, it will be passed to get_c_matrix() and the c_kwargs can be passed to this function. If sparse, the shape has to be (N_k_c, N_k) and it will be used directly. N_ker is the length of the effective kernel and N_k_c is the length of the spectrum (f_k) convolved. If None, the kernel is set to 1, i.e., do not do any convolution.

wave_gridarray-like, required

The grid on which f(lambda) will be projected, shape (N_k).

mask_trace_profilelist or ndarray

A list or array of 2-D boolean arrays of the pixel that need to be used for extraction, for each order on the detector. It has to have the same shape (N_ord, N, M) as trace_profile.

global_maskndarray, optional

Boolean mask of the detector pixels to mask for every extraction, e.g., bad pixels. Should not be related to a specific order (if so, use mask_trace_profile instead). Has shape (N, M).

orderslist, optional

List of orders considered. Default is [1, 2].

thresholdfloat, optional

The contribution of any order on a pixel is considered significant if its estimated spatial profile is greater than this threshold value. If it is not properly modeled (not covered by the wavelength grid), it will be masked. Default is 1e-3.

References

Darveau-Bernier 2022, PASP, DOI:10.1088/1538-3873/ac8a77

Attributes Summary

dtype
tikho_mat	Return the Tikhonov matrix, computing it if needed.

Methods Summary

__call__(data, error[, tikhonov, factor])	Extract underlying flux on the detector.
best_tikho_factor(tests, fit_mode)	Compute the best scale factor for Tikhonov regularization.
build_sys(data, error)	Build linear system arising from the logL maximisation.
compute_likelihood(spectrum, data, error)	Return the log likelihood associated with a particular spectrum.
compute_weights()	Compute integration weights.
estimate_tikho_factors(flux_estimate)	Estimate an initial guess of the Tikhonov factor.
get_attributes(*args[, i_order])	Return list of attributes.
get_detector_model(data, error)	Get the linear model of the detector pixel, B.dot(flux) = pixels.
get_mask_wave(i_order)	Generate mask bounded by limits of wavelength grid.
get_pixel_mapping(i_order[, error, quick])	Calculate the pixel mapping.
get_tikho_test_structure(data, error)	Structure to test different factors for Tikhonov regularization.
get_tikho_tests(tikho, factors)	Test different factors for Tikhonov regularization.
get_w(i_order)	Compute integration weights 'k' for each grid point and pixel 'i'.
grid_from_map([i_order])	Return the wavelength grid and the columns for a given order index.
precompute_detector_model(data, error, tikfac)	Return the matrices and mask needed to solve for the spectrum.
rebuild(spectrum[, fill_value])	Build current model image of the detector.
update_throughput(throughput)	Update internal throughput values.
wave_grid_c(i_order)	Return wave_grid for a given order constrained according to the i_bounds of that order.

Attributes Documentation

dtype = 'float64'

tikho_mat

Return the Tikhonov matrix, computing it if needed.

Methods Documentation

__call__(data, error, tikhonov=False, factor=None)

Extract underlying flux on the detector.

Performs an overlapping extraction of the form:

(B_T * B) * f = (data/sig)_T * B

where B is a matrix and f is an array. The matrix multiplication B * f is the 2D model of the detector. We want to solve for the array f. The elements of f are labelled by ‘k’. The pixels are labeled by ‘i’. Every pixel ‘i’ is covered by a set of ‘k’ for each order of diffraction.

Note

To be quicker, only specify the psf (p_list) in kwargs. There will be only one matrix multiplication:

(P/sig).(w.T.lambda.c_n)

Parameters:

data(N, M) array-like

A 2-D array of real values representing the detector image.

error(N, M) array-like

Estimate of the error on each pixel. Same shape as data.

tikhonovbool, optional

Whether to use Tikhonov extraction Default is False.

factorfloat, optional

The Tikhonov factor to use if tikhonov is True.

Returns:

spectrumarray[float]

Solution of the linear system (f_k).

best_tikho_factor(tests, fit_mode)

Compute the best scale factor for Tikhonov regularization.

The scale factor is determined by taking the factor giving the lowest reduced chi2 on the detector, the highest curvature of the l-curve or when the improvement on the chi2 (so the derivative of the chi2, ‘d_chi2’) reaches a certain threshold.

Parameters:

testsdict

Results of Tikhonov extraction tests for different factors. Must have the keys “factors” and “-logl”.

fit_modestr

Which mode is used to find the best Tikhonov factor. Options are ‘all’, ‘curvature’, ‘chi2’, or ‘d_chi2’. If ‘all’ is chosen, the best of the three other options will be selected.

Returns:

best_facfloat

The best Tikhonov factor.

build_sys(data, error)

Build linear system arising from the logL maximisation.

Parameters:

data(N, M) array-like

A 2-D array of real values representing the detector image.

error(N, M) array-like

Estimate of the error on each pixel.

Returns:

scipy.sparse.csr_array, array[float]

A, b from Ax = b being the system to solve.

compute_likelihood(spectrum, data, error)

Return the log likelihood associated with a particular spectrum.

Parameters:

spectrumarray[float] or callable

Flux as a function of wavelength if callable or array of flux values corresponding to self.wave_grid.

data(N, M) array-like

A 2-D array of real values representing the detector image.

error(N, M) array-like

Estimate of the error on each pixel. Same shape as data.

Returns:

array[float]

The log-likelihood of the spectrum.

compute_weights()

Compute integration weights.

The weights depend on the integration method used to solve the integral of the flux over a pixel and are encoded in the method get_w().

Returns:

weights, weights_k_idxlist

Lists of weights and corresponding grid indices

estimate_tikho_factors(flux_estimate)

Estimate an initial guess of the Tikhonov factor.

The output factor will be used to find the best Tikhonov factor. The flux_estimate is used to generate a factor_guess. The user should construct a grid with this output in log space, e.g., np.logspace(np.log10(flux_estimate)-4, np.log10(flux_estimate)+4, 9).

Parameters:

flux_estimatecallable

Estimate of the underlying flux (the solution f_k). Must be function of wavelengths and it will be projected on self.wave_grid.

Returns:

float

Estimated Tikhonov factor.

get_attributes(*args, i_order=None)

Return list of attributes.

Parameters:

*argsstr or list of str

All attributes to return.

i_orderNone or int, optional

Index of order to extract. If specified, it will be applied to all attributes in args, so it cannot be mixed with non-order dependent attributes.

Returns:

list

Result of [getattr(arg) for arg in args], with i_order indexing if provided.

get_detector_model(data, error)

Get the linear model of the detector pixel, B.dot(flux) = pixels.

Parameters:

dataarray-like

A 2-D array of real values representing the detector image, shape (N, M).

errorarray-like

Estimate of the error on each pixel, shape (N, M).

Returns:

B, pix_arrayarray[float]

From the linear equation:

B.dot(flux) = pix_array

get_mask_wave(i_order)

Generate mask bounded by limits of wavelength grid.

Parameters:

i_orderint

Order to select the wave_map on which a mask will be generated

Returns:

array[bool]

A mask with True where wave_map is outside the bounds of wave_grid

get_pixel_mapping(i_order, error=None, quick=False)

Calculate the pixel mapping.

Compute the matrix b_n = (P/sig).w.T.lambda.c_n, where P is the spatial profile matrix (diag), w is the integrations weights matrix, T is the throughput matrix (diag), lambda is the convolved wavelength grid matrix (diag), c_n is the convolution kernel. The model of the detector at order n (model_n) is given by the system: model_n = b_n.c_n.f, where f is the incoming flux projected on the wavelength grid. This method updates the b_n_list attribute.

Parameters:

i_orderint

Label of the order (depending on the initiation of the object).

errorarray-like or None, optional

Estimate of the error on each pixel. Same shape (N, M) as data. If None, the error is set to 1, which means the method will return b_n instead of b_n/sigma. Default is None.

quickbool, optional

If True, only perform one matrix multiplication instead of the whole system:

(P/sig).(w.T.lambda.c_n)

Returns:

array[float]

Sparse matrix of b_n coefficients.

get_tikho_test_structure(data, error)

Structure to test different factors for Tikhonov regularization.

Parameters:

data(N, M) array-like

A 2-D array of real values representing the detector image.

error(N, M) array-like

Estimate of the error on each pixel. Same shape as data.

Returns:

tikhoTikhonov

Instance of class with matrices pre-computed. Suitable for calling get_tikho_tests().

get_tikho_tests(tikho, factors)

Test different factors for Tikhonov regularization.

Parameters:

tikhoTikhonov

Instance of class with matrices pre-computed.

factors1D list or array-like

Factors to be tested.

Returns:

testsdict

Dictionary of the test results

get_w(i_order)

Compute integration weights ‘k’ for each grid point and pixel ‘i’.

These depend on the type of interpolation used, i.e., the order n.

Parameters:

i_orderint

Order to set the value of n in output arrays.

Returns:

w_narray

2D array of weights at this specific order n. The shape is given by: (number of pixels, max number of wavelengths covered by a pixel).

k_narray

2D array of the wavelength grid indices corresponding to the weights. Same shape as w_n.

grid_from_map(i_order=0)

Return the wavelength grid and the columns for a given order index.

Parameters:

i_orderint, optional

Order index to get the wavelength grid for. Default is 0.

Returns:

wave_gridarray[float]

Wavelength grid for the given order index.

icolarray[float]

Column indices for the wavelength grid.

precompute_detector_model(data, error, tikfac)

Return the matrices and mask needed to solve for the spectrum.

Parameters:

data(N, M) array-like

A 2-D array of real values representing the detector image.

error(N, M) array-like

Estimate of the error on each pixel. Same shape as data.

tikfacfloat

The Tikhonov factor to use.

Returns:

design_matrix_inv(N, N) array

The inverse of the design matrix M in M*f_k=b*(y/err).

b_matrix(N, M) array

The matrix b in M*f_k=b*(y/err).

rebuild(spectrum, fill_value=0.0)

Build current model image of the detector.

Parameters:

spectrumcallable or array-like

Flux as a function of wavelength if callable or array of flux values corresponding to self.wave_grid.

fill_valuefloat, optional

Pixel value where the detector is masked. Default is 0.0.

Returns:

array[float]

The modeled detector image.

update_throughput(throughput)

Update internal throughput values.

Parameters:

throughputarray[float] or callable

Throughput values for each order, given either as an array or as a callable function with self.wave_grid as input.

wave_grid_c(i_order)

Return wave_grid for a given order constrained according to the i_bounds of that order.

Parameters:

i_orderint

Order to select the wave_grid for.

Returns:

array[float]

wave_grid for the given order.