gelsa.optimize#
Functions
|
Build a block-diagonal discrete second-derivative operator D. |
|
Chi2, rms and pixel count over the target sources' footprint only. |
|
Build the full sparse Jacobian mapping each galaxy's spectrum to all masked pixels. |
|
|
|
|
|
|
|
|
|
|
|
Flag and remove outlier pixels based on residual/sigma, and rebuild the mask list to match. |
|
Extract spectra for all galaxies by building/solving the joint linear model across frames. |
|
|
|
Weighted least-squares solve of A @ x ~= b, with Tikhonov regularization. |
|
|
|
Boolean mask over the flat data vector selecting target-footprint pixels. |
|
Compute trace(Qinv @ M) without ever forming the n x n product. |
|
|
|
Rebuild per-detector pixel masks to account for newly invalidated pixels. |
- gelsa.optimize.build_smoothness_operator(galaxy_list)[source]#
Build a block-diagonal discrete second-derivative operator D.
D has one block per galaxy, in the same order as galaxy_list. Each block acts on that galaxy’s SED coefficients (its slice of the flat parameter vector x, as produced by unpack_extractions) and computes the discrete second derivative – the local curvature – of the spectrum along the wavelength grid:
(D_i @ x_i)[r] = x_i[r] - 2*x_i[r+1] + x_i[r+2]
Penalizing ||D @ x||^2 (instead of ||x||^2, as plain ridge does) therefore penalizes wiggliness/roughness in each galaxy’s spectrum rather than penalizing its overall flux level: a constant offset or a linear trend has zero second derivative, so those are left unpenalized while noisy high-frequency structure is suppressed.
The blocks are stacked block-diagonally (galaxy 0’s rows/columns, then galaxy 1’s, etc.) so that D never compares coefficients that belong to two different galaxies – this mirrors how compute_jacobian builds the Jacobian by horizontally stacking one block of columns per galaxy, and how unpack_extractions later splits x back into per-galaxy chunks in that same order.
- Parameters:
galaxy_list (list) – Galaxies in the same order used to build the Jacobian / x (each contributes len(galaxy.wavelength) columns to D).
- Returns:
D – n_cols equals len(x) (sum of len(gal.wavelength) over all galaxies). n_rows is the sum over galaxies of max(len(gal.wavelength) - 2, 0) – galaxies with fewer than 3 wavelength samples have no interior point to take a second difference of, so they contribute zero rows (no smoothness constraint on them, but they still need their columns present so the block widths add up to len(x)).
- Return type:
scipy.sparse csr array, shape (n_rows, n_cols)
- gelsa.optimize.trace_product(Qinv, M)[source]#
Compute trace(Qinv @ M) without ever forming the n x n product.
trace(Qinv @ M) = sum_ij Qinv_ij * M_ji, so the matmul – O(n^3) flops plus a second dense n x n allocation – is unnecessary. An elementwise product summed is O(n^2) dense, or O(nnz) when M is sparse. With n_params in the thousands that is the difference between seconds and milliseconds, and it is paid once per lambda per masking iteration, so it is worth the two lines.
Both matrices this is called with (Q_unreg, and the target-row Gram matrix) are symmetric Gram matrices, so the transpose is a formality – but it costs nothing (.T on a csr array is an O(1) view) and keeps the identity exact for any M.
- gelsa.optimize.solve(A, b, b_var, lam=0.001, compute_var=False, reg_operator=None, reg_floor=0.001, dof_rows=None)[source]#
Weighted least-squares solve of A @ x ~= b, with Tikhonov regularization.
Solves the regularized normal equations
(A^T Cinv A + lam_eff * P) x = A^T Cinv b
- where Cinv = diag(1/b_var) and P is the regularization matrix:
plain ridge (shrink-to-zero), if reg_operator is None: P = I, so the penalty term is lam_eff * ||x||^2.
smoothness (curvature) penalty, if reg_operator (e.g. from build_smoothness_operator) is given: P = reg_operator.T @ reg_operator, so the penalty term is lam_eff * ||reg_operator @ x||^2 – a small ridge floor is added on top, see reg_floor.
- Parameters:
A (sparse array) – Design matrix (Jacobian), rows = data points, columns = model parameters.
b (ndarray) – Data vector.
b_var (ndarray) – Per-data-point variance.
lam (float, optional) – Regularization strength; the matrix is scaled by the typical size of A’s entries (see lam_eff below) so that a given lam has a comparable effect regardless of the overall flux/units scale of A. lam <= 0 disables regularization.
compute_var (bool, optional) – If True, also return (an approximation of) the covariance of x, in the old post-hoc-smoothing style. See smooth above.
reg_operator (sparse array, optional) – Linear operator D such that ||D @ x||^2 is the penalty to minimize, instead of the plain ridge penalty ||x||^2.
reg_floor (float, optional) – Only used when reg_operator is given. Fraction (relative to the plain ridge strength) of a plain ||x||^2 penalty that is always added on top of the smoothness penalty. A pure second-difference penalty has a null space – a constant offset or a linear trend has zero curvature – so for a galaxy that’s only weakly constrained by data, the normal equations can come out singular without this floor to fall back on.
dof_rows (ndarray of bool, optional) – Boolean mask over the rows of A (i.e. over the entries of b) selecting a subset of the data. When given, the effective degrees of freedom are also computed restricted to those rows, and returned as dof_eff_rows. This is the diagonal sum of the hat matrix over that subset only – how much freedom the fit spends on those particular pixels – which is what pairs with a chi2 computed on the same subset to form a valid GCV-style score. The full parameter vector is still fit to all rows; only the accounting is restricted.
- Returns:
x (ndarray) – Solved parameter vector
estimator_var (ndarray)
chi2_reg (float) – The regularization penalty incurred by x, i.e. lam_eff * ||x||^2 (plain ridge) or lam_eff * ||reg_operator @ x||^2 (smoothness). Add this to the data-misfit chi2 to get the total penalized objective value.
dof_eff (float) – Effective degrees of freedom used by the fit (trace of the hat matrix). Use this in place of the nominal parameter count len(x) when forming a reduced chi2 for comparing different regularization strengths – see the comment where it’s computed below for why.
dof_eff_rows (float or None) – The same quantity restricted to the rows selected by dof_rows, or None when dof_rows was not given.
- gelsa.optimize.compute_jacobian(galaxy_list, frame_list, mask_list, frame_ordering, dx=1000, n_jobs=4, progress_hook=None)[source]#
Build the full sparse Jacobian mapping each galaxy’s spectrum to all masked pixels.
Computes the total number of unmasked pixels across all frames and detectors from mask_list, then builds a per-galaxy Jacobian block in parallel (via compute_jacobian_1) whose columns are the galaxy’s wavelength samples and whose rows are the concatenated masked pixels (ordered per frame_ordering). The per-galaxy blocks are horizontally stacked into a single sparse matrix so that a flat coefficient vector over all galaxies’ wavelengths maps directly onto the flat pixel vector produced by get_data_vec.
- Parameters:
galaxy_list (list or single galaxy) – Galaxy sources to build Jacobian columns for; normalized to a list via utils.ensurelist. Each contributes len(galaxy.wavelength) columns.
frame_list (list) – Frame/detector metadata used to evaluate the model image for each wavelength sample (passed through to compute_jacobian_1 / compute_jacobian_11 / evaluate_model).
mask_list (list) – Per-frame dict of {detector: (npix, mask)}; determines the total row count and which pixels are included.
frame_ordering (list) – Ordered (frame_i, det) pairs defining how per-detector pixel blocks are concatenated into rows.
dx (int, optional) – Sampling density controlling how densely each galaxy’s image is sampled per wavelength bin when evaluating the model (default 1000).
n_jobs (int, optional) – Number of parallel jobs used to compute each galaxy’s Jacobian block (default 4).
progress_hook (callable, optional) – Called with 1 each time a galaxy’s Jacobian block is done.
- Returns:
jac – Horizontally stacked sparse Jacobian of shape (total masked pixels, sum of galaxies’ wavelength lengths).
- Return type:
sparse array
- gelsa.optimize.update_masklist(valid, mask_list, frame_ordering)[source]#
Rebuild per-detector pixel masks to account for newly invalidated pixels.
Walks frame_ordering to slice the flat valid boolean array into per-(frame, detector) chunks matching each mask’s pixel count. For each detector, pixels that were already masked (mask <= -1) or are now invalid are excluded, and the surviving pixels are renumbered with consecutive indices (
0..npix_-1) so downstream code can index into the trimmed data/jac/var arrays. Each entry in mask_list is replaced with the updated(npix_, mask_)pair.- Parameters:
valid (ndarray of bool) – Flat validity flags, in the same concatenated order as frame_ordering, marking which pixels survive.
mask_list (list) – Per-frame dict/list of {det: (npix, mask)}, where mask maps original pixel positions to indices (-1 for already masked).
frame_ordering (list) – Ordered (frame_i, det) pairs describing how the flat valid array maps onto mask_list entries.
- Returns:
mask_list – The same structure with each (npix, mask) entry replaced by
(npix_, mask_)reflecting the newly masked/renumbered pixels.- Return type:
- gelsa.optimize.mask_jacobian(jac, data, var, resid, mask_list, frame_ordering, threshold=10)[source]#
Flag and remove outlier pixels based on residual/sigma, and rebuild the mask list to match.
Pixels with zero variance are treated as already masked. Among the remaining pixels, any whose residual exceeds threshold times the local sigma (sqrt(var)) are marked invalid. mask_list is then updated via update_masklist so its per-frame/detector indexing stays consistent with the newly masked pixels, and jac, data, and var are trimmed to the surviving (valid) rows.
- Parameters:
jac (ndarray) – Jacobian matrix, rows aligned with data/var/resid.
data (ndarray) – Data vector.
var (ndarray) – Per-pixel variance; entries <= 0 are treated as already masked.
resid (ndarray) – Residuals (data - model) used to compute the outlier statistic.
mask_list (list) – Per-frame mask info, as consumed/produced by update_masklist.
frame_ordering (list) – Ordering of (frame_i, det) pairs describing how resid maps onto mask_list entries.
threshold (float, optional) – Number of sigma above which a pixel is masked as an outlier (default 10).
- Returns:
jac, data, var (ndarray) – Inputs trimmed to the valid (non-masked) rows.
mask_list (list) – Updated mask list reflecting the newly masked pixels.
- gelsa.optimize.target_row_mask(frame_ordering, mask_list, mask_list_target, var=None)[source]#
Boolean mask over the flat data vector selecting target-footprint pixels.
Walks frame_ordering in the same order get_data_vec used to concatenate the data, and for each (frame, detector) block marks the surviving pixels that also fall inside the target sources’ footprint according to mask_list_target. Detectors absent from mask_list_target (no target lands on them) contribute an all-False block. If var is given, pixels with non-positive variance are dropped as well, matching what compute_chi2_target counts.
Factored out so that the target chi2 and the target-restricted effective dof are guaranteed to refer to the identical set of pixels – a GCV score built from a chi2 and a trace computed over different pixel sets would be silently wrong.
- Parameters:
frame_ordering (list) – Ordered (frame_i, det) pairs, as built by measure_everything.
mask_list (list) – Current per-frame {det: (npix, mask)} masks, i.e. the ones that match the current (post-sigma-clipping) data vector.
mask_list_target (list) – Per-frame {det: (npix, mask)} masks built from the target galaxies alone.
var (ndarray, optional) – Variance vector aligned with the data vector; when given, entries <= 0 are excluded.
- Returns:
sel – Length equals the current data vector length (sum of npix over frame_ordering).
- Return type:
ndarray of bool
- gelsa.optimize.compute_chi2_target(frame_ordering, mask_list, mask_list_target, data, model, var, sel=None)[source]#
Chi2, rms and pixel count over the target sources’ footprint only.
Pass sel (from target_row_mask) to reuse a mask already computed by the caller; otherwise it is derived here.
- gelsa.optimize.measure_everything(images, var_images, frame_list, target_galaxies, contaminant_galaxies=[], pixmask_list=None, jac=None, mask_list=None, dx=10000, n_jobs=4, masking_iterations=1, sigma_clip=5, smooth=True, reg_floor=0.001, reg_lambda_list=[1], compute_var=False, progress_hook=None)[source]#
Extract spectra for all galaxies by building/solving the joint linear model across frames.
Builds per-detector masks (via make_mask_1) if mask_list is not given, reconciles images and mask_list so both cover exactly the same frame/detector entries, and derives frame_ordering describing how per-detector pixels are concatenated. Builds the sparse Jacobian (via compute_jacobian) if not supplied, and assembles the flat data/variance vectors (via get_data_vec).
For each regularization strength in reg_lambda_list, iteratively solves the linear system (via solve), computing residuals and optionally re-masking outlier pixels (via mask_jacobian) for up to masking_iterations rounds before finalizing chi2/dof and unpacking the fitted coefficients into per-galaxy SED extractions.
- Parameters:
images (list) – Per-frame dict of {detector: image array} to extract from.
var_images (list) – Per-frame dict of {detector: variance array}, aligned with images.
frame_list (list) – Frame/detector metadata used to build masks and the Jacobian.
galaxies (list) – Galaxy sources to extract; each carries its own wavelength grid.
pixmask_list (list, optional) – Per-frame input pixel masks passed to make_mask_1 when mask_list is not supplied.
jac (sparse array, optional) – Precomputed Jacobian; if None, it is built via compute_jacobian.
reg_lambda_list (list, optional) – Regularization strengths to solve for; one result is returned per entry (default [1]).
smooth (bool, optional) – If True, regularize with a curvature (second-difference) penalty via build_smoothness_operator/solve’s reg_operator, instead of the plain ridge (shrink-to-zero) penalty (default False).
mask_list (list, optional) – Precomputed per-frame/detector masks; if None, built via make_mask_1 from pixmask_list.
dx (int, optional) – Sampling density passed to compute_jacobian (default 10000).
n_jobs (int, optional) – Number of parallel jobs used for mask building and Jacobian computation (default 4).
progress_hook (callable, optional) – Called with 1 as each parallel task finishes: one per frame for each mask pass (two passes when mask_list is None, one otherwise), then one per galaxy for the Jacobian when jac is None.
masking_iterations (int, optional) – Number of solve/re-mask rounds performed via mask_jacobian before the final solve for each lambda (default 1).
- Returns:
results – For each lambda in reg_lambda_list: galaxies_out is a copy of galaxies with .sed set to the extracted spectrum, and fit_info is a dict with frame_ordering, data, model, resid, var, chi2, dof, jac, mask_list, and lam.
No lambda is selected here – every entry in reg_lambda_list is returned and the choice is left to the caller. fit_info carries two generalized cross-validation scores for that purpose, each to be minimized over the returned results:
gcv_targ – computed over the target sources’ footprint only, and the one to select on. It scores the fit where it matters, on the pixels the extracted target spectrum is actually read from.
gcv – the same score over every fitted pixel. Kept as a diagnostic; it is dominated by contaminant and background pixels, so its minimum can sit at a different lambda.
See where they are computed below for the assumptions behind both.
- Return type:
list of (galaxies_out, fit_info)