"""
PCA utility functions for systematic noise removal in light curve residuals.
This module provides Principal Component Analysis (PCA) based functions
for processing light curve residuals and removing systematic components
in spectroscopic data analysis.
"""
import numpy as np
from sklearn.decomposition import PCA
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
[docs]
def linear_solver(
topca_l: np.ndarray,
res_l: np.ndarray
):
"""
Solve a linear system using Singular Value Decomposition (SVD).
This function solves the linear system Ax = b using SVD decomposition,
where A is the design matrix (res_l) and b is the target vector (topca_l).
The solution is computed in a numerically stable way using SVD.
Parameters:
topca_l (np.ndarray): Target vector b in the system Ax = b
res_l (np.ndarray): Design matrix A in the system Ax = b
Returns:
Results: The reconstructed array after solving the linear system
Error: Whether the linear system was solved successfully
"""
error = False
try:
# Perform SVD decomposition
uf, s, vh = np.linalg.svd(res_l, full_matrices=False)
# Compute intermediate matrices
uT = np.transpose(uf)
nD = np.diag(1 / s)
vHcT = np.transpose(vh.conj())
# Solve the system in steps:
# 1. U diag(s) Vh x = b <=> diag(s) Vh x = U.T b = c
c = np.dot(uT, topca_l)
# 2. diag(s) Vh x = c <=> Vh x = diag(1/s) c = w
w = np.dot(nD, c)
# 3. Vh x = w <=> x = Vh.H w
x = np.dot(vHcT, w)
# Reconstruct the array
return np.dot(res_l, x), error
except np.linalg.LinAlgError as e:
print(f"SVD computation failed: {str(e)}")
error = True
return None, error
[docs]
def trpca(
lcr_mask: np.ndarray,
lcrm_mask_nomask: np.ndarray,
good_pixel: np.ndarray,
ncomp: int,
smooth_on: bool,
smooth_size: int,
apply_standardize: bool,
model_reprocessing_type: str,
subtract_mean_spectrum: bool = False,
pca_mode: str = "spatial"
):
"""
Process light curve residuals by removing systematic components using PCA.
This function processes light curve residuals by:
1. Normalizing the data by subtracting column-wise means
2. Optionally applying PCA-based component removal
3. Updating the output array for good pixels
4. Optionally applying smoothing
Parameters:
lcr_mask (np.ndarray): 2D array of light curve residuals with telluric mask
lcrm_mask_nomask (np.ndarray): 2D array to store processed residuals
good_pixel (np.ndarray): Boolean mask or indices of good pixels
ncomp (int): Number of PCA components to use
smooth_on (bool): Whether to apply smoothing
smooth_size (int): Window size for smoothing filter
apply_standardize (bool): Whether to apply standardization
model_reprocessing_type (str): Processing mode, either "hard" or "soft"
subtract_mean_spectrum (bool): Whether to subtract mean spectrum
pca_mode (str): PCA mode - "spatial" (samples=images, features=pixels) or
"temporal" (samples=pixels, features=images). Default: "spatial"
Returns:
np.ndarray: Processed light curve residuals
error: Whether the processing was successful
"""
error = False
# Normalize the light curve residuals
lcrm = np.copy(lcr_mask)
means = np.mean(lcrm, axis=0)
lcrm = lcrm - np.expand_dims(means, axis=0)
# Center data to column
if subtract_mean_spectrum:
mean_spectrum = np.mean(lcrm, axis=1)
lcrm = lcrm - np.expand_dims(mean_spectrum, axis=1)
# Apply PCA-based component removal if in "hard" mode
if model_reprocessing_type == "hard":
try:
# Prepare data based on PCA mode
# Spatial mode (default): samples=images, features=pixels -> matrix shape (n_pixels, n_images)
# Temporal mode: samples=pixels, features=images -> matrix shape (n_images, n_pixels)
if pca_mode == "temporal":
# Transpose for temporal PCA: now rows are images, columns are pixels
pca_input_matrix = lcrm.T
else:
# Spatial PCA (default): rows are pixels, columns are images
pca_input_matrix = lcrm
# Create pipeline with optional standardization
if apply_standardize:
transformer = Pipeline([
('scaler', StandardScaler()),
('pca', PCA(n_components=ncomp, svd_solver="full"))
])
else:
transformer = Pipeline([
('pca', PCA(n_components=ncomp, svd_solver="full"))
])
comps = transformer.fit_transform(pca_input_matrix)
ones_column = np.ones((comps.shape[0], 1))
comps_with_ones = np.hstack((comps, ones_column))
except Exception as e: # except (np.linalg.LinAlgError, ValueError) as e:
print(f"Error applying PCA-based component removal: {str(e)}")
error = True
return None, error
#
array_reconstruct, error_trpca = linear_solver(pca_input_matrix, comps_with_ones)
if error_trpca:
error = True
print(f"PCA processing failed in solver")
return None, error
#
diff_pca = pca_input_matrix - array_reconstruct
# Transform back to original orientation if temporal mode was used
if pca_mode == "temporal":
lcrm_pca = diff_pca.T
else:
lcrm_pca = diff_pca
else:
lcrm_pca = lcrm
# Update output array for good pixels
lcrm_mask_nomask[good_pixel, :] = lcrm_pca
# Apply smoothing if requested
if smooth_on:
# Note: Smoothing implementation is currently disabled
# This could be reimplemented using scipy.ndimage.uniform_filter
pass
return lcrm_pca, error