mccd.mccd

MCCD CLASS.

This module contains the main MCCD class.

Authors

Tobias Liaudat <tobias.liaudat@cea.fr>

Jerome Bonnin <https://github.com/jerome-bonnin>

Morgan Schmitz <https://github.com/MorganSchmitz>

mccd_quickload(path)[source]

Load pre-fitted MCCD model.

(saved with mccd.quicksave()).

Parameters

path (str) – Path to the npy file containing the saved MCCD model.

Returns

loaded_model – The MCCD model.

Return type

MCCD class

class MCCD(n_comp_loc, d_comp_glob, d_hyb_loc=2, min_d_comp_glob=None, upfact=1, ksig_loc=1.0, ksig_glob=1.0, rmse_thresh=1.25, ccd_star_thresh=0.15, n_scales=3, ksig_init=1.0, filters=None, verbose=2, fp_geometry='CFIS')[source]

Bases: object

Multi-CCD Resolved Components Analysis.

Parameters

  • n_comp_loc (int) – Number of components to learn for each CCD. The interpretation may depend on the loc_model parameter of the fit() function.

  • d_comp_glob (int) – Degree of polynomial components for the global model.

  • d_hyb_loc (int) – Degree of the polynomial component for the local part in case the local model used is hybrid. Default is 2.

  • min_d_comp_glob (int or None) – The minimum degree of the polynomial for the global component. For example, if the paramter is set to 1, the polynomials of degree 0 and 1 will be excluded from the global polynomial variations. None means that we are not excluding any degree. Default is None.

  • upfact (int) – Upsampling factor. Default is 1 (no superresolution).

  • ksig_loc (float) – Value of \(K_{\sigma}^{Loc}\) for the thresholding in wavelet (Starlet by default) domain (taken to be \(K\sigma\), where \(\sigma\) is the estimated noise standard deviation.). It is used for the thresholding of the local eigenPSF \(S_{K}\) update. Default is 1.

  • ksig_glob (float) – Value of \(K_{\sigma}^{Glob}\) for the thresholding in Starlet domain (taken to be \(k\sigma\), where \(\sigma\) is the estimated noise standard deviation.). It is used for the thresholding of the global eigenPSF \(\tilde{S}\) update. Default is 1.

  • rmse_thresh (float) – Parameter concerning the CCD outlier rejection. Once the PSF model is calculated we perform an outlier check on the training stars. We divide each star in two parts with a given circle. The inner part corresponds to the most of the PSF/star energy while the outer part corresponds to the observation background. The outer part is used to calculate the noise level and the inner part to calculate the model residual (star observation - PSF model reconstruction). If the RMSE error of the residual divided by the noise level is over the rmse_thresh the star will be considered an outlier. A perfect reconstruction would have rmse_thresh equal to 1. Default is 1.25.

  • ccd_star_thresh (float) – Parameter concerning the CCD outlier rejection. If the percentage of outlier stars in a single CCD is bigger than ccd_star_thresh, the CCD is considered to be an outlier. In this case, the CCD is rejected from the PSF model. A value lower than 0 means that no outlier rejection will be done. Default is 0.15.

  • n_scales (int) – Number of wavelet (Default Starlet) scales to use for the sparsity constraint. Default is 3. Unused if filters are provided.

  • ksig_init (float) – Similar to ksig, for use when estimating shifts and noise levels, as it might be desirable to have it set higher than ksig. Unused if shifts are provided when running RCA.fit(). Default is 5.

  • filters (numpy.ndarray) – Optional filters to the transform domain wherein eigenPSFs are assumed to be sparse; convolution by them should amount to applying \(\Phi\). Optional; if not provided, the Starlet transform with n_scales scales will be used.

  • verbose (bool or int) – If True, will only output RCA-specific lines to stdout. If verbose is set to 2, will run ModOpt’s optimization algorithms in verbose mode. Default to 2.

  • fp_geometry (str) – Geometry of the focal plane. It defines the transformations to use to go from local to global coordinates. Default is CFIS. Available options are CFIS, EUCLID.

quicksave(path)[source]

Save fitted model.

Save fitted MCCD model for later use. Ideally, you would probably want to store the whole MCCD instance, though this might mean storing a lot of data you are not likely to use if you do not alter the fit that was already performed. Stored models can be loaded with mccd.mccd_quickload().

Parameters

path (str) – Path to where the fitted MCCDF model should be saved.

fit(obs_data, obs_pos, ccd_list, obs_weights=None, SNR_weight_list=None, S=None, VT=None, Pi=None, alpha=None, shifts=None, sigs=None, psf_size=6, psf_size_type='R2', flux=None, nb_iter=1, nb_iter_glob=2, nb_iter_loc=2, nb_subiter_S_loc=100, nb_reweight=0, nb_subiter_A_loc=500, nb_subiter_S_glob=30, nb_subiter_A_glob=200, n_eigenvects=5, loc_model='hybrid', pi_degree=2, graph_kwargs={})[source]

Fits MCCD to observed star field.

Parameters

  • obs_data (list of numpy.ndarray) – Observed data (each element of the list being one CCD).

  • obs_pos (list of numpy.ndarray) – Corresponding positions (global coordinate system).

  • ccd_list (list of numpy.ndarray) – List containing the ccd_ids of each set of observations, positions and weights. It is of utmost importance that the ccd_list contains the ccd_id in the same order as in the other lists. Ex: obs_data[0] is the data from the ccd ccd_list[0].

  • obs_weights (list of numpy.ndarray) – Corresponding weights. Can be either one per observed star, or contain pixel-wise values. Masks can be handled via binary weights. Default is None (in which case no weights are applied). Note if fluxes and shifts are not provided, weights will be ignored for their estimation. Noise level estimation only removes bad pixels (with weight strictly equal to 0) and otherwise ignores weights.

  • SNR_weight_list (list of floats) – List of weights to be used on each star. They can probably be based on a specific function of the SNR and should represent the degree of significance we give to a specific star. The values should be around 1. Default is None meaning that no weights will be used.

  • S (list of numpy.ndarray) – First guess (or warm start) eigenPSFs \(S\) (last matrix is global). Default is None.

  • VT (list of numpy.ndarray) – Matrices of concatenated eigenvectors of the different graph Laplacians. Default is None.

  • Pi (list of numpy.ndarray) – Matrices of polynomials in positions. Default is None.

  • alpha (list numpy.ndarray) – First guess (or warm start) weights \(\alpha\), after factorization by VT (last matrix is global). Default is None.

  • shifts (list of numpy.ndarray) – Corresponding sub-pixel shifts. Default is None; will be estimated from observed data if not provided.

  • sigs (numpy.ndarray) – Estimated noise levels. Default is None; will be estimated from data if not provided.

  • psf_size (float) – Approximate expected PSF size in psf_size_type; will be used for the size of the Gaussian window for centroid estimation. psf_size_type determines the convention used for this size (default is FWHM). Ignored if shifts are provided. Default is 'R2' of 6.

  • psf_size_type (str) – Can be any of 'R2', 'fwhm' or 'sigma', for the size defined from quadrupole moments, full width at half maximum (e.g. from SExtractor) or 1-sigma width of the best matching 2D Gaussian. 'sigma' is the value outputed from Galsim’s HSM adaptive moment estimator. Default is 'R2'.

  • flux (list of numpy.ndarray) – Flux levels. Default is None; will be estimated from data if not provided.

  • nb_iter (int) – Number of overall iterations (i.e. of alternations). Note the weights and global components do not get updated the last time around, so they actually get nb_iter-1 updates. Paramter \(l_{max}\) on the MCCD article pseudo-algorithm. Default is 1.

  • nb_iter_glob (int) – Number of iterations on the global model estimation. The times we go trough the global \(S,A\) updates. Paramter \(n_G\) on the MCCD article pseudo-algorithm. Default value is 2.

  • nb_iter_loc (int) – Number of iterations on the local model estimation. The times we go trough the local \(S,A\) updates. Paramter \(n_L\) on the MCCD article pseudo-algorithm. Default value is 2.

  • nb_subiter_S_loc (int) – Number of iterations when solving the optimization problem (III) concerning the local eigenPSFs \(S_{k}\). Default is 100.

  • nb_subiter_A_loc (int) – Number of iterations when solving the optimization problem (IV) concerning the local weights \(\alpha_{k}\). Default is 500.

  • nb_subiter_S_glob (int) – Number of iterations when solving the optimization problem (I) concerning the global eigenPSFs \(\tilde{S}\). Default is 30.

  • nb_subiter_A_glob (int) – Number of iterations when solving the optimization problem (II) concerning the global weights \(\tilde{\alpha}\). Default is 200.

  • nb_reweight (int) – Number of reweightings to apply during \(S\) updates. See equation (33) in RCA paper. Default is 0.

  • n_eigenvects (int) – Maximum number of eigenvectors to consider per \((e,a)\) couple. Default is None; if not provided, all eigenvectors will be considered, which can lead to a poor selection of graphs, especially when data is undersampled. Ignored if VT and alpha are provided.

  • loc_model (str) – Defines the type of local model to use, it can be: 'rca', 'poly' or 'hybrid'. Thus defining the MCCD-RCA, MCCD-POL and MCCD-HYB. When MCCD-POL is used, n_comp_loc should be used as the d_comp_glob (max degree of the polynomial) for the local model. When MCCD-HYB is used, n_comp_loc should be used as in MCCD-RCA, the number of graph-based eigenPSFs. The max local polynomial degree is set to 2.

  • pi_degree (int) – Maximum degree of polynomials in Pi. Default is 2. Ignored if Pi is provided.

  • graph_kwargs (dict) – List of optional kwargs to be passed on to the utils.GraphBuilder().

remove_ccd_from_model(ccd_idx)[source]

Remove ccd from the trained model.

remove_outlier_ccds()[source]

Remove all CCDs with outliers.

Reminder: the outlier rejection is done on the train stars. We will reject a CCD if the percentage of outlier stars in a single CCD is bigger than ccd_star_thresh.

They outlier threshold is based in rmse_thresh. A perfect reconstruction would have rmse_thresh equal to 1.

interpolate_psf_pipeline(test_pos, ccd_n, centroid=None)[source]

Estimate PSF at desired position with the required centroid.

This function is a consequence of following the requirements needed to use the MCCD model in a pipeline. (ShapePipe now but could be adapted)

The returned PSFs should be normalized with unitary flux and with the required centroid.

Please note the differences between the pixel conventions between this package and Galsim. In the first one the position of a pixel is in its lower left corner and the indexing starts at zero (numpy style). Galsim’s convention is that the position is in the center of the pixel and the indexation starts at one.

Returns the model in “regular” format, (n_stars, n_pixels, n_pixels).

Parameters

  • test_pos (numpy.ndarray) – Positions where the PSF should be estimated. Should be in the same format (coordinate system, units, etc.) as the obs_pos fed to MCCD.fit().

  • ccd_n (int) – ccd_id of the positions to be tested.

  • centroid (list of floats, [xc, yc]) – Required centroid for the output PSF. It has to be specified following the MCCD pixel convention. Default will be the vignet centroid.

Returns

PSFs – Returns the interpolated PSFs in regular format.

Return type

numpy.ndarray

estimate_psf(test_pos, ccd_n, n_loc_neighbors=15, n_glob_neighbors=15, rbf_function='thin_plate', apply_degradation=False, shifts=None, flux=None, sigmas=None, upfact=None, mccd_debug=False, global_pol_interp=None)[source]

Estimate and return PSF at desired positions.

Returns the model in “regular” format, (n_stars, n_pixels, n_pixels).

Parameters

  • test_pos (numpy.ndarray) – Positions where the PSF should be estimated. Should be in the same format (coordinate system, units, etc.) as the obs_pos fed to MCCD.fit().

  • ccd_n (int) – ccd_id of the positions to be tested.

  • n_loc_neighbors (int) – Number of neighbors for the local model to use for RBF interpolation. Default is 15.

  • n_glob_neighbors (int) – Number of neighbors for the global model to use for RBF interpolation. Default is 15.

  • rbf_function (str) – Type of RBF kernel to use. Default is 'thin_plate'. Check scipy RBF documentation for more models.

  • apply_degradation (bool) – Whether PSF model should be degraded (shifted and resampled on coarse grid), for instance for comparison with stars. If True, expects shifts to be provided. Needed if pixel validation against stars will be required. Default is False.

  • shifts (numpy.ndarray) – Intra-pixel shifts to apply if apply_degradation is set to True. Needed to match the observed stars in case of pixel validation.

  • flux (numpy.ndarray) – Flux levels by which reconstructed PSF will be multiplied if provided. For pixel validation with stars if apply_degradation is set to True.

  • sigmas (numpy.ndarray) – Sigmas (shapes) are used to have a better estimate of the size of the lanczos interpolant that will be used when performing the intra-pixel shift. Default is None, where a predefined value is used for the interpolant size.

  • upfact (int) – Upsampling factor; default is None, in which case that of the MCCD instance will be used.

  • mccd_debug (bool) – Debug option. It returns the model plus the local and the global reconstruction components. Default is False.

  • global_pol_interp (numpy.ndarray) – If is None, the global interpolation is done with th RBF interpolation as in the local model. If is not None, the global interpolation is done directly using position polynomials. In this case, global_pol_interp should be the normalized Pi interpolation matrix. Default is None.

Returns

PSFs – Returns the interpolated PSFs in regular format.

Return type

numpy.ndarray

validation_stars(test_stars, test_pos, test_masks=None, ccd_id=None, mccd_debug=False, response_flag=False, global_pol_interp=None)[source]

Match PSF model to stars.

The match is done in flux, shift and pixel sampling - for validation tests. Returns the matched PSFs’ stamps. Returned matched PSFs will be None if there is no model on that specific CCD due to the lack of training stars.

Parameters

  • test_stars (numpy.ndarray) – Star stamps to be used for comparison with the PSF model. Should be in “rca” format, i.e. with axises (n_pixels, n_pixels, n_stars).

  • test_pos (numpy.ndarray) – Their corresponding positions in global coordinate system.

  • test_masks (numpy.ndarray) – Masks to be used on the star stamps. If None, all the pixels will be used. Default is None.

  • ccd_id (int) – The corresponding ccd_id (ie ccd number corresponding to the megacam geometry). Do not mistake for ccd_idx (index).

  • mccd_debug (bool) – Debug option. It returns the local and the global reconstruction components.

  • response_flag (bool) – Response option. True if in response mode.

  • global_pol_interp (Position pols of numpy.ndarray or None) – If is None, the global interpolation is done with th RBF interpolation as in the local model. If is not None, the global interpolation is done directly using position polynomials. In this case, it should be the normalized Pi interpolation matrix.