diff --git a/shapepipe/modules/mccd_fit_runner.py b/shapepipe/modules/mccd_fit_runner.py index cdca57948..403294a61 100644 --- a/shapepipe/modules/mccd_fit_runner.py +++ b/shapepipe/modules/mccd_fit_runner.py @@ -9,20 +9,32 @@ """ from shapepipe.modules.module_decorator import module_runner -from shapepipe.modules.MCCD_package import shapepipe_auxiliary_mccd as aux_mccd +from shapepipe.modules.mccd_package import shapepipe_auxiliary_mccd\ + as aux_mccd import mccd -@module_runner(input_module=['mccd_preprocessing_runner'], version='1.0', - file_pattern=['train_star_selection'], - file_ext=['.fits'], numbering_scheme='-0000000', - depends=['numpy', 'mccd', 'galsim'], run_method='parallel') -def mccd_fit_runner(input_file_list, run_dirs, file_number_string, - config, w_log): +@module_runner( + input_module=['mccd_preprocessing_runner'], + version='1.0', + file_pattern=['train_star_selection'], + file_ext=['.fits'], + numbering_scheme='-0000000', + depends=['numpy', 'mccd', 'galsim'], + run_method='parallel' +) +def mccd_fit_runner( + input_file_list, + run_dirs, + file_number_string, + config, + module_config_sec, + w_log +): # Recover the MCCD config file and its params - config_file_path = config.getexpanded('MCCD', 'CONFIG_PATH') - mccd_mode = config.get('MCCD', 'MODE') - verbose = config.getboolean('MCCD', 'VERBOSE') + config_file_path = config.getexpanded(module_config_sec, 'CONFIG_PATH') + mccd_mode = config.get(module_config_sec, 'MODE') + verbose = config.getboolean(module_config_sec, 'VERBOSE') # Parse MCCD config file mccd_parser = mccd.auxiliary_fun.MCCDParamsParser(config_file_path) @@ -35,11 +47,18 @@ def mccd_fit_runner(input_file_list, run_dirs, file_number_string, if mccd_mode == 'FIT': aux_mccd.mccd_fit_pipeline( - trainstar_path, file_number_string, mccd_parser, output_dir, - verbose, saving_name) + trainstar_path=trainstar_path, + file_number_string=file_number_string, + mccd_parser=mccd_parser, + output_dir=output_dir, + verbose=verbose, + saving_name=saving_name, + w_log=w_log + ) else: - raise ValueError('''mccd_fit_runner should be called when the - MODE is "FIT".''') + raise ValueError( + "mccd_fit_runner should be called when the MODE is 'FIT'." + ) return None, None diff --git a/shapepipe/modules/mccd_fit_val_runner.py b/shapepipe/modules/mccd_fit_val_runner.py index 0391edea1..189b0fd06 100644 --- a/shapepipe/modules/mccd_fit_val_runner.py +++ b/shapepipe/modules/mccd_fit_val_runner.py @@ -9,20 +9,31 @@ """ from shapepipe.modules.module_decorator import module_runner -from shapepipe.modules.MCCD_package import shapepipe_auxiliary_mccd as aux_mccd +from shapepipe.modules.mccd_package import shapepipe_auxiliary_mccd as aux_mccd import mccd -@module_runner(input_module=['mccd_preprocessing_runner'], version='1.0', - file_pattern=['train_star_selection', 'test_star_selection'], - file_ext=['.fits', '.fits'], numbering_scheme='-0000000', - depends=['numpy', 'mccd', 'galsim'], run_method='parallel') -def mccd_fit_val_runner(input_file_list, run_dirs, file_number_string, - config, w_log): +@module_runner( + input_module=['mccd_preprocessing_runner'], + version='1.0', + file_pattern=['train_star_selection', 'test_star_selection'], + file_ext=['.fits', '.fits'], + numbering_scheme='-0000000', + depends=['numpy', 'mccd', 'galsim'], + run_method='parallel' +) +def mccd_fit_val_runner( + input_file_list, + run_dirs, + file_number_string, + config, + module_config_sec, + w_log +): # Recover the MCCD config file and its params - config_file_path = config.getexpanded('MCCD', 'CONFIG_PATH') - mccd_mode = config.get('MCCD', 'MODE') - verbose = config.getboolean('MCCD', 'VERBOSE') + config_file_path = config.getexpanded(module_config_sec, 'CONFIG_PATH') + mccd_mode = config.get(module_config_sec, 'MODE') + verbose = config.getboolean(module_config_sec, 'VERBOSE') # Parse MCCD config file mccd_parser = mccd.auxiliary_fun.MCCDParamsParser(config_file_path) @@ -45,7 +56,9 @@ def mccd_fit_val_runner(input_file_list, run_dirs, file_number_string, mccd_parser=mccd_parser, output_dir=output_dir, verbose=verbose, - saving_name=fit_saving_name) + saving_name=fit_saving_name, + w_log=w_log + ) # Fitted model is found in the output directory mccd_model_path = output_dir + fit_saving_name + file_number_string \ @@ -58,10 +71,13 @@ def mccd_fit_val_runner(input_file_list, run_dirs, file_number_string, output_dir=output_dir, file_number_string=file_number_string, w_log=w_log, - val_saving_name=val_saving_name) + val_saving_name=val_saving_name + ) else: - raise ValueError('''mccd_fit_val_runner should be called when the - MODE is "FIT_VALIDATION".''') + raise ValueError( + "mccd_fit_val_runner should be called when the MODE" + + " is 'FIT_VALIDATION'." + ) return None, None diff --git a/shapepipe/modules/mccd_interp_runner.py b/shapepipe/modules/mccd_interp_runner.py index 6b5d7ad96..f4b42398e 100644 --- a/shapepipe/modules/mccd_interp_runner.py +++ b/shapepipe/modules/mccd_interp_runner.py @@ -9,34 +9,49 @@ """ from shapepipe.modules.module_decorator import module_runner -from shapepipe.modules.MCCD_package import mccd_interpolation_script\ +from shapepipe.modules.mccd_package import mccd_interpolation_script\ as mccd_interp -from shapepipe.modules.MCCD_package import shapepipe_auxiliary_mccd as aux_mccd +from shapepipe.modules.mccd_package import shapepipe_auxiliary_mccd\ + as aux_mccd import os -@module_runner(input_module=['setools_runner'], - version='1.0', - file_pattern=['galaxy_selection'], - file_ext=['.npy'], - depends=['numpy', 'astropy', 'galsim', 'sqlitedict', 'mccd'], - run_method='parallel') -def mccd_interp_runner(input_file_list, run_dirs, file_number_string, - config, w_log): - - mode = config.getexpanded('MCCD_INTERP_RUNNER', 'MODE') - pos_params = config.getlist('MCCD_INTERP_RUNNER', 'POSITION_PARAMS') - get_shapes = config.getboolean('MCCD_INTERP_RUNNER', 'GET_SHAPES') +@module_runner( + input_module=['setools_runner'], + version='1.0', + file_pattern=['galaxy_selection'], + file_ext=['.npy'], + depends=['numpy', 'astropy', 'galsim', 'sqlitedict', 'mccd'], + run_method='parallel' +) +def mccd_interp_runner( + input_file_list, + run_dirs, + file_number_string, + config, + module_config_sec, + w_log +): + + mode = config.getexpanded(module_config_sec, 'MODE') + pos_params = config.getlist(module_config_sec, 'POSITION_PARAMS') + get_shapes = config.getboolean(module_config_sec, 'GET_SHAPES') mccd_model_extension = '.npy' output_dir = run_dirs['output'] if mode == 'CLASSIC': - psf_model_dir = config.getexpanded('MCCD_INTERP_RUNNER', - 'PSF_MODEL_DIR') - psf_model_pattern = config.get('MCCD_INTERP_RUNNER', - 'PSF_MODEL_PATTERN') - psf_separator = config.get('MCCD_INTERP_RUNNER', - 'PSF_MODEL_SEPARATOR') + psf_model_dir = config.getexpanded( + module_config_sec, + 'PSF_MODEL_DIR' + ) + psf_model_pattern = config.get( + module_config_sec, + 'PSF_MODEL_PATTERN' + ) + psf_separator = config.get( + module_config_sec, + 'PSF_MODEL_SEPARATOR' + ) # psfcat_path, galcat_path = input_file_list galcat_path = input_file_list[0] @@ -52,8 +67,10 @@ def mccd_interp_runner(input_file_list, run_dirs, file_number_string, if not os.path.exists(psf_model_path): error_msg = "The corresponding PSF model was not not found." - w_log.info("Error message: {}.\n On catalog with id: {}.".format( - error_msg, file_number_string)) + w_log.info( + f"Error message: {error_msg}.\n On catalog with" + + f" id: {file_number_string}." + ) return None, None @@ -64,37 +81,50 @@ def mccd_interp_runner(input_file_list, run_dirs, file_number_string, pos_params=pos_params, ccd_id=int(ccd_id), saving_path=saving_path, - get_shapes=get_shapes) + get_shapes=get_shapes + ) if output_msg is not None: - w_log.info("Error message: {}.\n On catalog with id: {}.".format( - output_msg, file_number_string)) + w_log.info( + f"Error message: {output_msg}.\n On catalog with" + + f" id: {file_number_string}." + ) elif mode == 'MULTI-EPOCH': - psf_model_dir = config.getexpanded('MCCD_INTERP_RUNNER', - 'PSF_MODEL_DIR') - psf_model_pattern = config.get('MCCD_INTERP_RUNNER', - 'PSF_MODEL_PATTERN') - f_wcs_path = config.getexpanded('MCCD_INTERP_RUNNER', - 'ME_LOG_WCS') + psf_model_dir = config.getexpanded( + module_config_sec, + 'PSF_MODEL_DIR' + ) + psf_model_pattern = config.get( + module_config_sec, + 'PSF_MODEL_PATTERN' + ) + f_wcs_path = config.getexpanded( + module_config_sec, + 'ME_LOG_WCS' + ) galcat_path = input_file_list[0] - inst = mccd_interp.MCCDinterpolator(None, - galcat_path, - output_dir, - file_number_string, - w_log, - pos_params, - get_shapes) + inst = mccd_interp.MCCDinterpolator( + None, + galcat_path, + output_dir, + file_number_string, + w_log, + pos_params, + get_shapes + ) inst.process_me(psf_model_dir, psf_model_pattern, f_wcs_path) elif mode == 'VALIDATION': - ValueError('''MODE has to be in MULTI-EPOCH or CLASSIC. - For validation use MCCD validation runner''') + ValueError( + "MODE has to be in MULTI-EPOCH or CLASSIC. For validation" + + " use MCCD validation runner." + ) else: - ValueError('MODE has to be in : [CLASSIC, MULTI-EPOCH]') + ValueError("MODE has to be in : [CLASSIC, MULTI-EPOCH]") return None, None diff --git a/shapepipe/modules/mccd_merge_starcat_runner.py b/shapepipe/modules/mccd_merge_starcat_runner.py index 36318e3ff..5819776ec 100644 --- a/shapepipe/modules/mccd_merge_starcat_runner.py +++ b/shapepipe/modules/mccd_merge_starcat_runner.py @@ -15,13 +15,23 @@ from shapepipe.modules.module_decorator import module_runner -@module_runner(input_module=['mccd_fit_val_runner'], version='1.0', - file_pattern=['validation_psf'], - file_ext=['.fits'], numbering_scheme='-0000000', - depends=['numpy', 'astropy'], - run_method='serial') -def mccd_merge_starcat_runner(input_file_list, run_dirs, file_number_string, - config, w_log): +@module_runner( + input_module=['mccd_fit_val_runner'], + version='1.0', + file_pattern=['validation_psf'], + file_ext=['.fits'], + numbering_scheme='-0000000', + depends=['numpy', 'astropy'], + run_method='serial' +) +def mccd_merge_starcat_runner( + input_file_list, + run_dirs, + file_number_string, + config, + module_config_sec, + w_log +): w_log.info('Merging validation results..') hdu_table = 1 @@ -69,16 +79,19 @@ def mccd_merge_starcat_runner(input_file_list, run_dirs, file_number_string, # Pixel mse calculation pix_val = np.sum((stars - psfs) ** 2) pix_sum = np.sum((stars - psfs)) - masked_diffs = np.array([(_star - _psf)[my_mask] for _star, _psf in - zip(stars, psfs)]) + masked_diffs = np.array( + [(_star - _psf)[my_mask] for _star, _psf in zip(stars, psfs)] + ) masked_pix_val = np.sum(masked_diffs ** 2) masked_pix_sum = np.sum(masked_diffs) # Star noise variance (using masked stars) star_noise_var_val = np.array( - [np.var(_star[~my_mask]) for _star in stars]) - res_var_val = np.array([np.var(_star - _psf) for _star, _psf in - zip(stars, psfs)]) + [np.var(_star[~my_mask]) for _star in stars] + ) + res_var_val = np.array( + [np.var(_star - _psf) for _star, _psf in zip(stars, psfs)] + ) # Variance of the model # (Residual variance - Star variance (using masked stars)) @@ -90,16 +103,21 @@ def mccd_merge_starcat_runner(input_file_list, run_dirs, file_number_string, stars_norm_vals = np.sqrt(np.sum(stars ** 2, axis=(1, 2))) psfs_norm_vals = np.sqrt(np.sum(psfs ** 2, axis=(1, 2))) # Select non zero stars & psfs - non_zero_elems = np.logical_and((psfs_norm_vals != 0), - (stars_norm_vals != 0)) + non_zero_elems = np.logical_and( + (psfs_norm_vals != 0), + (stars_norm_vals != 0) + ) # Calculate the filtered mse calculation - pix_filt_val = np.sum((stars[non_zero_elems] - - psfs[non_zero_elems]) ** 2) + pix_filt_val = np.sum( + (stars[non_zero_elems] - psfs[non_zero_elems]) ** 2 + ) # Calculate the normalized (& filtered) mse calculation stars_norm_vals = stars_norm_vals[non_zero_elems].reshape(-1, 1, 1) psfs_norm_vals = psfs_norm_vals[non_zero_elems].reshape(-1, 1, 1) - pix_norm_val = np.sum((stars[non_zero_elems] / stars_norm_vals - - psfs[non_zero_elems] / psfs_norm_vals) ** 2) + pix_norm_val = np.sum( + (stars[non_zero_elems] / stars_norm_vals - + psfs[non_zero_elems] / psfs_norm_vals) ** 2 + ) # Calculate sizes filt_size = stars[non_zero_elems].size regular_size = stars.size @@ -132,12 +150,16 @@ def mccd_merge_starcat_runner(input_file_list, run_dirs, file_number_string, ra += list(starcat_j[hdu_table].data['RA_LIST'][:]) dec += list(starcat_j[hdu_table].data['DEC_LIST'][:]) except Exception: - ra += list(np.zeros( - starcat_j[hdu_table].data['GLOB_POSITION_IMG_LIST'][:, 0].shape, - dtype=int)) - dec += list(np.zeros( - starcat_j[hdu_table].data['GLOB_POSITION_IMG_LIST'][:, 0].shape, - dtype=int)) + ra += list( + np.zeros(starcat_j[hdu_table].data[ + 'GLOB_POSITION_IMG_LIST'][:, 0].shape, + dtype=int) + ) + dec += list( + np.zeros(starcat_j[hdu_table].data[ + 'GLOB_POSITION_IMG_LIST'][:, 0].shape, + dtype=int) + ) # shapes (convert sigmas to R^2) g1_psf += list(starcat_j[hdu_table].data['PSF_MOM_LIST'][:, 0]) @@ -160,56 +182,63 @@ def mccd_merge_starcat_runner(input_file_list, run_dirs, file_number_string, # Stats about the merging def rmse_calc(values, sizes): return np.sqrt( - np.nansum(np.array(values)) / np.nansum(np.array(sizes))) + np.nansum(np.array(values)) / np.nansum(np.array(sizes)) + ) def rmse_calc_2(values, sizes): return np.sqrt( - np.nansum(np.array(values) ** 2) / np.nansum(np.array(sizes))) + np.nansum(np.array(values) ** 2) / np.nansum(np.array(sizes)) + ) def mean_calc(values, sizes): - return np.nansum(np.array(values)) / np.nansum( - np.array(sizes)) + return np.nansum( + np.array(values)) / np.nansum(np.array(sizes) + ) def std_calc(values): return np.nanstd(np.array(values)) # Regular pixel RMSE tot_pixel_rmse = rmse_calc(pixel_mse, size_mse) - # print('Regular Total pixel RMSE = %.5e\n' % tot_pixel_rmse) - w_log.info('MCCD_merge_starcat: Regular Total pixel RMSE = %.5e\n' % ( - tot_pixel_rmse)) + w_log.info( + f"MCCD_merge_starcat: Regular Total pixel RMSE =" + + f" {tot_pixel_rmse:.5e}\n" + ) # Regular Total pixel mean tot_pixel_mean = mean_calc(pixel_sum, size_mse) - # print('Regular Total pixel mean = %.5e\n' % tot_pixel_mean) - w_log.info('MCCD_merge_starcat: Regular Total pixel mean = %.5e\n' % ( - tot_pixel_mean)) + w_log.info( + f"MCCD_merge_starcat: Regular Total pixel mean =" + + f" {tot_pixel_mean:.5e}\n" + ) # Regular Total MASKED pixel RMSE tot_masked_pixel_rmse = rmse_calc(masked_pixel_mse, masked_size) - # print('Regular Total MASKED pixel RMSE = %.5e\n' % tot_masked_pixel_rmse) w_log.info( - 'MCCD_merge_starcat: Regular Total MASKED pixel RMSE = %.5e\n' % ( - tot_masked_pixel_rmse)) + f"MCCD_merge_starcat: Regular Total MASKED pixel RMSE =" + + f" {tot_masked_pixel_rmse:.5e}\n" + ) # Regular Total MASKED pixel mean tot_masked_pixel_mean = mean_calc(masked_pixel_sum, masked_size) - # print('Regular Total MASKED pixel mean = %.5e\n' % tot_masked_pixel_mean) w_log.info( - 'MCCD_merge_starcat: Regular Total MASKED pixel mean = %.5e\n' % ( - tot_masked_pixel_mean)) + f"MCCD_merge_starcat: Regular Total MASKED pixel mean =" + + f" {tot_masked_pixel_mean:.5e}\n" + ) # Normalized pixel RMSE tot_pix_norm_rmse = rmse_calc(pix_norm_mse, size_norm_mse) - # print('Normalized Total pixel RMSE = %.5e\n' % tot_pix_norm_rmse) - w_log.info('MCCD_merge_starcat: Normalized Total pixel RMSE = %.5e\n' % ( - tot_pix_norm_rmse)) + w_log.info( + "MCCD_merge_starcat: Normalized Total pixel RMSE =" + + f" {tot_pix_norm_rmse:.5e}\n" + ) # Normalized filtered pixel RMSE tot_pix_filt_rmse = rmse_calc(pix_filt_mse, size_filt_mse) - # print('Filtered Total pixel RMSE = %.5e\n' % tot_pix_filt_rmse) - w_log.info('MCCD_merge_starcat: Filtered Total pixel RMSE = %.5e\n' % ( - tot_pix_filt_rmse)) + w_log.info( + "MCCD_merge_starcat: Filtered Total pixel RMSE =" + + f" {tot_pix_filt_rmse:.5e}\n" + ) concat_model_var = np.concatenate(np.array(model_var)) concat_star_noise_var = np.concatenate(np.array(star_noise_var)) @@ -218,23 +247,21 @@ def std_calc(values): mean_model_var = mean_calc(concat_model_var, model_var_size) std_model_var = std_calc(concat_model_var) rmse_model_var = rmse_calc_2(concat_model_var, model_var_size) - # print("MCCD-RCA variance:\n Mean Variance= %.5e\n " - # "Std Variance= %.5e\n RMSE Variance= %.5e\n" % ( - # mean_model_var, std_model_var, rmse_model_var)) - w_log.info("MCCD-RCA variance:\n Mean Variance= %.5e\n " - "Std Variance= %.5e\n RMSE Variance= %.5e\n" % - (mean_model_var, std_model_var, rmse_model_var)) + w_log.info( + f"MCCD-RCA variance:\n Mean Variance= {mean_model_var:.5e}\n" + + f"Std Variance= {std_model_var:.5e}\n" + + f"RMSE Variance= {rmse_model_var:.5e}\n" + ) # Star Noise Variance mean_star_var = mean_calc(concat_star_noise_var, star_noise_size) std_star_var = std_calc(concat_star_noise_var) rmse_star_var = rmse_calc_2(concat_star_noise_var, star_noise_size) - # print("Masked stars variance:\n Mean Variance= %.5e\n" - # "Std Variance= %.5e\n RMSE Variance= %.5e\n" % ( - # mean_star_var, std_star_var, rmse_star_var)) - w_log.info("Masked stars variance:\n Mean Variance= %.5e\n" - "Std Variance= %.5e\n RMSE Variance= %.5e\n" % - (mean_star_var, std_star_var, rmse_star_var)) + w_log.info( + f"Masked stars variance:\n Mean Variance= {mean_star_var:.5e}\n" + + f"Std Variance= {std_star_var:.5e}\n" + + f"RMSE Variance= {rmse_star_var:.5e}\n" + ) # Mask and transform to numpy arrays flagmask = np.abs(np.array(flag_star) - 1) * np.abs(np.array(flag_psf) - 1) @@ -246,38 +273,48 @@ def std_calc(values): star_r2 = np.array(size)[flagmask.astype(bool)] rmse, mean, std_dev = stats_calculator(star_e1, psf_e1) - # print("Moment residual e1:\n Mean= %.5e\n " - # "Std Dev= %.5e\n RMSE= %.5e\n" % (mean, std_dev, rmse)) - w_log.info("Moment residual e1:\n Mean= %.5e\n" - "Std Dev= %.5e\n RMSE= %.5e\n" % (mean, std_dev, rmse)) + w_log.info( + f"Moment residual e1:\n Mean= {mean:.5e}\nStd Dev= {std_dev:.5e}\n" + + f" RMSE= {rmse:.5e}\n" + ) rmse, mean, std_dev = stats_calculator(star_e2, psf_e2) - # print("Moment residual e2:\n Mean= %.5e\n" - # "Std Dev= %.5e\n RMSE= %.5e\n" % (mean, std_dev, rmse)) - w_log.info("Moment residual e2:\n Mean= %.5e\n" - "Std Dev= %.5e\n RMSE= %.5e\n" % (mean, std_dev, rmse)) + w_log.info( + f"Moment residual e2:\n Mean= {mean:.5e}\nStd Dev= {std_dev:.5e}\n" + + f"RMSE= {rmse:.5e}\n" + ) rmse, mean, std_dev = stats_calculator(star_r2, psf_r2) - # print("Moment residual R2:\n Mean= %.5e\n" - # "Std Dev= %.5e\n RMSE= %.5e\n" % (mean, std_dev, rmse)) - w_log.info("Moment residual R2:\n Mean= %.5e\n" - "Std Dev= %.5e\n RMSE= %.5e\n" % (mean, std_dev, rmse)) + w_log.info( + f"Moment residual R2:\n Mean= {mean:.5e}\nStd Dev= {std_dev:.5e}\n" + + f"RMSE= {rmse:.5e}\n" + ) - print('MCCD: Number of stars: %d' % (star_e1.shape[0])) - w_log.info('MCCD: Number of stars: %d' % (star_e1.shape[0])) + print(f"MCCD: Number of stars: {star_e1.shape[0]:d}") + w_log.info(f"MCCD: Number of stars: {star_e1.shape[0]:d}") - output = sc.FITSCatalog(output_dir + '/full_starcat-0000000.fits', - open_mode=sc.BaseCatalog.OpenMode.ReadWrite, - SEx_catalog=True) + output = sc.FITSCatalog( + output_dir + '/full_starcat-0000000.fits', + open_mode=sc.BaseCatalog.OpenMode.ReadWrite, + SEx_catalog=True + ) # convert back to sigma for consistency - data = {'X': x, 'Y': y, 'RA': ra, 'DEC': dec, - 'E1_PSF_HSM': g1_psf, 'E2_PSF_HSM': g2_psf, - 'SIGMA_PSF_HSM': np.sqrt(size_psf), - 'E1_STAR_HSM': g1, 'E2_STAR_HSM': g2, - 'SIGMA_STAR_HSM': np.sqrt(size), - 'FLAG_PSF_HSM': flag_psf, 'FLAG_STAR_HSM': flag_star, - 'CCD_NB': ccd_nb} + data = { + 'X': x, + 'Y': y, + 'RA': ra, + 'DEC': dec, + 'E1_PSF_HSM': g1_psf, + 'E2_PSF_HSM': g2_psf, + 'SIGMA_PSF_HSM': np.sqrt(size_psf), + 'E1_STAR_HSM': g1, + 'E2_STAR_HSM': g2, + 'SIGMA_STAR_HSM': np.sqrt(size), + 'FLAG_PSF_HSM': flag_psf, + 'FLAG_STAR_HSM': flag_star, + 'CCD_NB': ccd_nb + } print('Writing full catalog...') output.save_as_fits(data, sex_cat_path=input_file_list[0][0]) diff --git a/shapepipe/modules/MCCD_package/__init__.py b/shapepipe/modules/mccd_package/__init__.py similarity index 85% rename from shapepipe/modules/MCCD_package/__init__.py rename to shapepipe/modules/mccd_package/__init__.py index 180f073f7..afcdeb80a 100644 --- a/shapepipe/modules/MCCD_package/__init__.py +++ b/shapepipe/modules/mccd_package/__init__.py @@ -1,6 +1,6 @@ # -*- coding: utf-8 -*- -"""MCCD auxiliary functions package. +"""MCCD AUXILIARY FUNCTIONS PACKAGE. This module contains auxiliary functions for the MCCD package. diff --git a/shapepipe/modules/MCCD_package/mccd_interpolation_script.py b/shapepipe/modules/mccd_package/mccd_interpolation_script.py similarity index 76% rename from shapepipe/modules/MCCD_package/mccd_interpolation_script.py rename to shapepipe/modules/mccd_package/mccd_interpolation_script.py index 2c3e2816f..b173d3679 100644 --- a/shapepipe/modules/MCCD_package/mccd_interpolation_script.py +++ b/shapepipe/modules/mccd_package/mccd_interpolation_script.py @@ -59,8 +59,10 @@ def interp_MCCD(mccd_model_path, positions, ccd): # Positions to global coordinates loc2glob = mccd.mccd_utils.Loc2Glob() - glob_pos = np.array([loc2glob.loc2glob_img_coord(_ccd, _pos[0], _pos[1]) - for _ccd, _pos in zip(ccd_list, positions)]) + glob_pos = np.array([ + loc2glob.loc2glob_img_coord(_ccd, _pos[0], _pos[1]) + for _ccd, _pos in zip(ccd_list, positions) + ]) # Interpolate the model PSFs = mccd_instance.estimate_psf(glob_pos, ccd) @@ -102,8 +104,16 @@ class MCCDinterpolator(object): """ - def __init__(self, dotpsf_path, galcat_path, output_path, img_number, - w_log, pos_params=None, get_shapes=True): + def __init__( + self, + dotpsf_path, + galcat_path, + output_path, + img_number, + w_log, + pos_params=None, + get_shapes=True + ): # Path to PSFEx output file self._dotpsf_path = dotpsf_path # Path to catalog containing galaxy positions @@ -127,9 +137,10 @@ def __init__(self, dotpsf_path, galcat_path, output_path, img_number, # CosmoStat's ShapePipe) if pos_params: if not len(pos_params) == 2: - raise ValueError('{} position parameters were passed on; ' - 'there should be exactly two.' - ''.format(len(pos_params))) + raise ValueError( + f'{len(pos_params)} position parameters were passed on;' + + f' there should be exactly two.' + ) self._pos_params = pos_params else: self._pos_params = None @@ -143,8 +154,10 @@ def _get_position_parameters(self): """ dotpsf = sc.FITSCatalog(self._dotpsf_path) dotpsf.open() - self._pos_params = [dotpsf.get_header()['POLNAME1'], - dotpsf.get_header()['POLNAME2']] + self._pos_params = [ + dotpsf.get_header()['POLNAME1'], + dotpsf.get_header()['POLNAME2'] + ] dotpsf.close() def _get_galaxy_positions(self): @@ -159,15 +172,19 @@ def _get_galaxy_positions(self): self.gal_pos = np.array( [[x, y] for x, y in zip( galcat.get_data()[self._pos_params[0]], - galcat.get_data()[self._pos_params[1]])]) + galcat.get_data()[self._pos_params[1]] + )] + ) except KeyError as detail: # extract erroneous position parameter from original exception err_pos_param = detail.args[0][4:-15] - pos_param_err = ('Required position parameter ' + err_pos_param + - 'was not found in galaxy catalog. Leave ' - 'pos_params (or EXTRA_CODE_OPTION) blank to ' - 'read them from .psf file.') + pos_param_err = ( + 'Required position parameter ' + err_pos_param + + 'was not found in galaxy catalog. Leave ' + 'pos_params (or EXTRA_CODE_OPTION) blank to ' + 'read them from .psf file.' + ) raise KeyError(pos_param_err) galcat.close() @@ -177,28 +194,38 @@ def _get_psfshapes(self): if import_fail: raise ImportError('Galsim is required to get shapes information') - psf_moms = [hsm.FindAdaptiveMom(Image(psf), strict=False) - for psf in self.interp_PSFs] - - self.psf_shapes = np.array([[moms.observed_shape.g1, - moms.observed_shape.g2, - moms.moments_sigma, - int(bool(moms.error_message))] for moms in - psf_moms]) + psf_moms = [ + hsm.FindAdaptiveMom(Image(psf), strict=False) + for psf in self.interp_PSFs + ] + + self.psf_shapes = np.array([ + [ + moms.observed_shape.g1, + moms.observed_shape.g2, + moms.moments_sigma, + int(bool(moms.error_message)) + ] + for moms in psf_moms + ]) def _write_output(self): """ Save computed PSFs to fits file. """ - output = sc.FITSCatalog(self._output_path + self._img_number + '.fits', - open_mode=sc.BaseCatalog.OpenMode.ReadWrite, - SEx_catalog=True) + output = sc.FITSCatalog( + self._output_path + self._img_number + '.fits', + open_mode=sc.BaseCatalog.OpenMode.ReadWrite, + SEx_catalog=True + ) if self._compute_shape: - data = {'VIGNET': self.interp_PSFs, - 'E1_PSF_HSM': self.psf_shapes[:, 0], - 'E2_PSF_HSM': self.psf_shapes[:, 1], - 'SIGMA_PSF_HSM': self.psf_shapes[:, 2], - 'FLAG_PSF_HSM': self.psf_shapes[:, 3].astype(int)} + data = { + 'VIGNET': self.interp_PSFs, + 'E1_PSF_HSM': self.psf_shapes[:, 0], + 'E2_PSF_HSM': self.psf_shapes[:, 1], + 'SIGMA_PSF_HSM': self.psf_shapes[:, 2], + 'FLAG_PSF_HSM': self.psf_shapes[:, 3].astype(int) + } else: data = {'VIGNET': self.interp_PSFs} output.save_as_fits(data, sex_cat_path=self._galcat_path) @@ -216,15 +243,22 @@ def _get_starshapes(self, star_vign): masks = np.zeros_like(star_vign) masks[np.where(star_vign == -1e30)] = 1 - star_moms = [ - hsm.FindAdaptiveMom(Image(star), badpix=Image(mask), strict=False) - for star, mask in zip(star_vign, masks)] - - self.star_shapes = np.array([[moms.observed_shape.g1, - moms.observed_shape.g2, - moms.moments_sigma, - int(bool(moms.error_message))] for moms - in star_moms]) + star_moms = [hsm.FindAdaptiveMom( + Image(star), + badpix=Image(mask), + strict=False) + for star, mask in zip(star_vign, masks) + ] + + self.star_shapes = np.array([ + [ + moms.observed_shape.g1, + moms.observed_shape.g2, + moms.moments_sigma, + int(bool(moms.error_message)) + ] + for moms in star_moms + ]) def _get_psfexcatdict(self, psfex_cat_path): """ Get data from PSFEx .cat file. @@ -266,14 +300,16 @@ def _write_output_validation(self, star_dict, psfex_cat_dict): open_mode=sc.BaseCatalog.OpenMode.ReadWrite, SEx_catalog=True) - data = {'E1_PSF_HSM': self.psf_shapes[:, 0], - 'E2_PSF_HSM': self.psf_shapes[:, 1], - 'SIGMA_PSF_HSM': self.psf_shapes[:, 2], - 'FLAG_PSF_HSM': self.psf_shapes[:, 3].astype(int), - 'E1_STAR_HSM': self.star_shapes[:, 0], - 'E2_STAR_HSM': self.star_shapes[:, 1], - 'SIGMA_STAR_HSM': self.star_shapes[:, 2], - 'FLAG_STAR_HSM': self.star_shapes[:, 3].astype(int)} + data = { + 'E1_PSF_HSM': self.psf_shapes[:, 0], + 'E2_PSF_HSM': self.psf_shapes[:, 1], + 'SIGMA_PSF_HSM': self.psf_shapes[:, 2], + 'FLAG_PSF_HSM': self.psf_shapes[:, 3].astype(int), + 'E1_STAR_HSM': self.star_shapes[:, 0], + 'E2_STAR_HSM': self.star_shapes[:, 1], + 'SIGMA_STAR_HSM': self.star_shapes[:, 2], + 'FLAG_STAR_HSM': self.star_shapes[:, 3].astype(int) + } data = {**data, **star_dict} data['ACCEPTED'] = np.ones_like(data['NUMBER'], dtype='int16') @@ -345,38 +381,53 @@ def _interpolate_me(self): # self._dot_psf_pattern + '-' + exp_name + '-' + str(ccd) +\ # '.psf' mccd_model_path = self._dot_psf_dir + '/' +\ - self._dot_psf_pattern + '-' + exp_name +\ - '.npy' + self._dot_psf_pattern + '-' + exp_name + '.npy' + ind_obj = np.where(cat.get_data(hdu_index)['CCD_N'] == ccd)[0] obj_id = all_id[ind_obj] gal_pos = np.array( self._f_wcs_file[exp_name][ccd]['WCS'].all_world2pix( self.gal_pos[:, 0][ind_obj], - self.gal_pos[:, 1][ind_obj], 0)).T + self.gal_pos[:, 1][ind_obj], + 0 + ) + ).T self.interp_PSFs = interp_MCCD(mccd_model_path, gal_pos, ccd) # self.interp_PSFs = interpsfex( # dot_psf_path, gal_pos, self._star_thresh, self._chi2_thresh) - if isinstance(self.interp_PSFs, - str) and self.interp_PSFs == NOT_ENOUGH_STARS: - self._w_log.info('Not enough stars find in the ccd' - ' {} of the exposure {}. Object inside' - ' this ccd will lose an' - 'epoch.'.format(ccd, exp_name)) + if ( + isinstance(self.interp_PSFs, str) + and (self.interp_PSFs == NOT_ENOUGH_STARS) + ): + self._w_log.info( + f'Not enough stars find in the ccd {ccd} of the' + + f' exposure {exp_name}. Object inside this ' + + f'ccd will lose an epoch.' + ) continue - if isinstance(self.interp_PSFs, - str) and self.interp_PSFs == BAD_CHI2: - self._w_log.info('Bad chi2 for the psf model in the ccd' - ' {} of the exposure {}. Object inside' - ' this ccd will lose an' - 'epoch.'.format(ccd, exp_name)) + + if ( + isinstance(self.interp_PSFs, str) + and (self.interp_PSFs == BAD_CHI2) + ): + self._w_log.info( + f'Bad chi2 for the psf model in the ccd {ccd} of the' + + f' exposure {exp_name}. Object inside this ccd' + + f' will lose an epoch.' + ) + continue - if isinstance(self.interp_PSFs, - str) and self.interp_PSFs == FILE_NOT_FOUND: - self._w_log.info('Psf model file {} not found. ' - 'Object inside this ccd will lose' - 'an epoch'.format(self._dotpsf_path)) + + if ( + isinstance(self.interp_PSFs, str) + and (self.interp_PSFs == FILE_NOT_FOUND) + ): + self._w_log.info( + f'Psf model file {self._dotpsf_path} not found.' + + f' Object inside this ccd will lose an epoch.' + ) continue if array_psf is None: diff --git a/shapepipe/modules/MCCD_package/mccd_plot_utilities.py b/shapepipe/modules/mccd_package/mccd_plot_utilities.py similarity index 75% rename from shapepipe/modules/MCCD_package/mccd_plot_utilities.py rename to shapepipe/modules/mccd_package/mccd_plot_utilities.py index a50920b90..f03f1b956 100644 --- a/shapepipe/modules/MCCD_package/mccd_plot_utilities.py +++ b/shapepipe/modules/mccd_package/mccd_plot_utilities.py @@ -84,8 +84,14 @@ def MegaCamFlip_2(xbins, ybins, ccd_nb, nb_pixel): return xbins, ybins -def MeanShapesPlot(ccd_maps, filename, title='', colorbar_ampl=1., wind=None, - cmap='bwr'): +def MeanShapesPlot( + ccd_maps, + filename, + title='', + colorbar_ampl=1., + wind=None, + cmap='bwr' +): r"""Plot meanshapes from ccd maps.""" # colorbar amplitude if wind is None: @@ -102,22 +108,33 @@ def MeanShapesPlot(ccd_maps, filename, title='', colorbar_ampl=1., wind=None, axes.flat[j].axis('off') for ccd_nb, ccd_map in enumerate(ccd_maps): ax = axes.flat[MegaCamPos(ccd_nb)] - im = ax.imshow(ccd_map.T, cmap=cmap, interpolation='Nearest', - vmin=vmin, vmax=vmax) + im = ax.imshow( + ccd_map.T, + cmap=cmap, + interpolation='Nearest', + vmin=vmin, + vmax=vmax + ) ax.set_xticks([]) ax.set_yticks([]) - ax.set_title('rmse=%.3e' % (np.sqrt(np.nanmean(ccd_map ** 2))), size=8) + ax.set_title(f'rmse={np.sqrt(np.nanmean(ccd_map ** 2)):.3e}', size=8) plt.suptitle(title, size=20) # TODO: fix title fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax) - plt.savefig('{}.png'.format(filename)) + plt.savefig(f'{filename}.png') plt.close() -def plot_meanshapes(starcat_path, output_path, nb_pixel, w_log, - remove_outliers=False, plot_meanshapes=True, - plot_histograms=True): +def plot_meanshapes( + starcat_path, + output_path, + nb_pixel, + w_log, + remove_outliers=False, + plot_meanshapes=True, + plot_histograms=True +): r"""Plot meanshapes and histograms.""" # READ FULL STARCAT starcat = fits.open(starcat_path, memmap=False) @@ -146,14 +163,27 @@ def plot_meanshapes(starcat_path, output_path, nb_pixel, w_log, all_X = starcat[2].data['X'] all_Y = starcat[2].data['Y'] + # Remove stars/PSFs where the measured size is zero + # Sometimes the HSM shape measurement gives objects with measured + # size equals to zero without an error Flag. + bad_stars = (abs(all_star_shapes[2, :]) < 0.1) + bad_psfs = (abs(all_psf_shapes[2, :]) < 0.1) + size_mask = np.abs(bad_stars) * np.abs(bad_psfs) + # Remove outlier stars/PSFs + all_star_shapes = all_star_shapes[:, ~size_mask] + all_psf_shapes = all_psf_shapes[:, ~size_mask] + all_CCDs = all_CCDs[~size_mask] + all_X = all_X[~size_mask] + all_Y = all_Y[~size_mask] + flagmask = flagmask[~size_mask] + if remove_outliers: shape_std_max = 5. # Outlier rejection based on the size R2_thresh = shape_std_max * np.std(all_star_shapes[2, :]) + np.mean( all_star_shapes[2, :]) bad_stars = (abs(all_star_shapes[2, :]) > R2_thresh) - bad_stars_idx = np.nonzero(bad_stars)[0] - w_log.info('Nb of outlier stars: %d' % (np.sum(bad_stars))) + w_log.info(f'Nb of outlier stars: {np.sum(bad_stars):d}') # Remove outlier PSFs all_star_shapes = all_star_shapes[:, ~bad_stars] all_psf_shapes = all_psf_shapes[:, ~bad_stars] @@ -162,12 +192,18 @@ def plot_meanshapes(starcat_path, output_path, nb_pixel, w_log, all_Y = all_Y[~bad_stars] flagmask = flagmask[~bad_stars] - w_log.info('TOTAL e1 residual RMSE: %.6e\n' % ( - np.sqrt(np.mean((all_star_shapes[0, :] - all_psf_shapes[0, :]) ** 2)))) - w_log.info('TOTAL e2 residual RMSE: %.6e\n' % ( - np.sqrt(np.mean((all_star_shapes[1, :] - all_psf_shapes[1, :]) ** 2)))) - w_log.info('TOTAL R2 residual RMSE: %.6e\n' % ( - np.sqrt(np.mean((all_star_shapes[2, :] - all_psf_shapes[2, :]) ** 2)))) + e1_res_rmse = np.sqrt( + np.mean((all_star_shapes[0, :] - all_psf_shapes[0, :]) ** 2) + ) + e2_res_rmse = np.sqrt( + np.mean((all_star_shapes[1, :] - all_psf_shapes[1, :]) ** 2) + ) + R2_res_rmse = np.sqrt( + np.mean((all_star_shapes[2, :] - all_psf_shapes[2, :]) ** 2) + ) + w_log.info(f"TOTAL e1 residual RMSE: {e1_res_rmse:.6e}\n") + w_log.info(f"TOTAL e2 residual RMSE: {e2_res_rmse:.6e}\n") + w_log.info(f"TOTAL R2 residual RMSE: {R2_res_rmse:.6e}\n") # CCDs x star/model x (e1,e2,R2,nstars) x xpos x ypos ccd_maps = np.ones((40, 2, 4) + nb_pixel) * np.nan @@ -210,85 +246,125 @@ def plot_meanshapes(starcat_path, output_path, nb_pixel, w_log, np.abs(np.nanmin(ccd_maps[:, :, 0]))) vmin = -vmax wind = [vmin, vmax] + title = ( + f"e_1 (stars), std={np.nanstd(ccd_maps[:, 0, 0]):.5e}\n" + + f"vmax={np.nanmax(abs(ccd_maps[:, 0, 0])):.4e}" + ) MeanShapesPlot( - ccd_maps[:, 0, 0], output_path + 'e1s', - "e_1 (stars), std=%.5e\nvmax=%.4e" % ( - np.nanstd(ccd_maps[:, 0, 0]), - np.nanmax(abs(ccd_maps[:, 0, 0]))), - wind=wind) - + ccd_maps[:, 0, 0], + output_path + 'e1s', + title, + wind=wind + ) + + title = ( + f"e_1 (model), std={np.nanstd(ccd_maps[:, 1, 0]):.5e}\n" + + f"vmax={np.nanmax(abs(ccd_maps[:, 1, 0])):.4e}" + ) MeanShapesPlot( - ccd_maps[:, 1, 0], output_path + 'e1m', - 'e_1 (model), std=%.5e\nvmax=%.4e' % - (np.nanstd(ccd_maps[:, 1, 0]), - np.nanmax(abs(ccd_maps[:, 1, 0]))), - wind=wind) + ccd_maps[:, 1, 0], + output_path + 'e1m', + title, + wind=wind + ) if auto_colorbar: wind = None e1_res = ccd_maps[:, 0, 0] - ccd_maps[:, 1, 0] e1_res = e1_res[~np.isnan(e1_res)] rmse_e1 = np.sqrt(np.mean(e1_res ** 2)) - w_log.info('Bins: e1 residual RMSE: %.6f\n' % rmse_e1) + w_log.info(f"Bins: e1 residual RMSE: {rmse_e1:.6f}\n") vmax = np.nanmax(abs(ccd_maps[:, 0, 0] - ccd_maps[:, 1, 0])) vmin = -vmax wind = [vmin, vmax] + title = ( + f"e_1 res, rmse={rmse_e1:.5e}\nvmax={vmax:.4e} , " + + f"std={np.nanstd(ccd_maps[:, 0, 0] - ccd_maps[:, 1, 0]):.5e}" + ) MeanShapesPlot( ccd_maps[:, 0, 0] - ccd_maps[:, 1, 0], output_path + 'e1res', - 'e_1 res, rmse=%.5e\nvmax=%.4e , std=%.5e' % ( - rmse_e1, vmax, - np.nanstd(ccd_maps[:, 0, 0] - ccd_maps[:, 1, 0])), + title, wind=wind, - colorbar_ampl=colorbar_ampl) + colorbar_ampl=colorbar_ampl + ) # e_2 - vmax = max(np.nanmax(ccd_maps[:, :, 1]), - np.abs(np.nanmin(ccd_maps[:, :, 1]))) + vmax = max( + np.nanmax(ccd_maps[:, :, 1]), + np.abs(np.nanmin(ccd_maps[:, :, 1])) + ) vmin = -vmax wind = [vmin, vmax] - MeanShapesPlot(ccd_maps[:, 0, 1], output_path + 'e2s', - 'e_2 (stars), std=%.5e\nvmax=%.4e' % - (np.nanstd(ccd_maps[:, 0, 1]), - np.nanmax(abs(ccd_maps[:, 0, 1]))), wind=wind) - MeanShapesPlot(ccd_maps[:, 1, 1], output_path + 'e2m', - 'e_2 (model), std=%.5e\nvmax=%.4e' % - (np.nanstd(ccd_maps[:, 1, 1]), - np.nanmax(abs(ccd_maps[:, 1, 1]))), wind=wind) + title = ( + f"e_2 (stars), std={np.nanstd(ccd_maps[:, 0, 1]):.5e}\n" + + f"vmax={np.nanmax(abs(ccd_maps[:, 0, 1])):.4e}" + ) + MeanShapesPlot( + ccd_maps[:, 0, 1], + output_path + 'e2s', + title, + wind=wind + ) + title = ( + f"e_2 (model), std={np.nanstd(ccd_maps[:, 1, 1]):.5e}\n" + + f"vmax={np.nanmax(abs(ccd_maps[:, 1, 1])):.4e}" + ) + MeanShapesPlot( + ccd_maps[:, 1, 1], + output_path + 'e2m', + title, + wind=wind + ) if auto_colorbar: wind = None colorbar_ampl = 1. + e2_res = ccd_maps[:, 0, 1] - ccd_maps[:, 1, 1] e2_res = e2_res[~np.isnan(e2_res)] rmse_e2 = np.sqrt(np.mean(e2_res ** 2)) - w_log.info('Bins: e2 residual RMSE: %.6f\n' % rmse_e2) + w_log.info(f"Bins: e2 residual RMSE: {rmse_e2:.6f}\n") vmax = np.nanmax(abs(ccd_maps[:, 0, 1] - ccd_maps[:, 1, 1])) vmin = -vmax wind = [vmin, vmax] + title = ( + f"e_2 res, rmse={rmse_e2:.5e}\nvmax={vmax:.4e} , " + + f"std={np.nanstd(ccd_maps[:, 0, 1] - ccd_maps[:, 1, 1]):.5e}" + ) MeanShapesPlot( - ccd_maps[:, 0, 1] - ccd_maps[:, 1, 1], output_path + 'e2res', - 'e_2 res, rmse=%.5e\nvmax=%.4e , std=%.5e' % ( - rmse_e2, vmax, - np.nanstd(ccd_maps[:, 0, 1] - ccd_maps[:, 1, 1])), - wind=wind, colorbar_ampl=colorbar_ampl) + ccd_maps[:, 0, 1] - ccd_maps[:, 1, 1], + output_path + 'e2res', + title, + wind=wind, + colorbar_ampl=colorbar_ampl + ) # R^2 wind = [0, np.nanmax(ccd_maps[:, :, 2])] colorbar_ampl = 1 + title = ( + f"R_2 (stars), std={np.nanstd(ccd_maps[:, 0, 2]):.5e}\n" + + f"vmax={np.nanmax(abs(ccd_maps[:, 0, 2])):.4e}" + ) MeanShapesPlot( - ccd_maps[:, 0, 2], output_path + 'R2s', - 'R_2 (stars), std=%.5e\nvmax=%.4e' % ( - np.nanstd(ccd_maps[:, 0, 2]), - np.nanmax(abs(ccd_maps[:, 0, 2]))), - wind=wind, cmap='Reds') - + ccd_maps[:, 0, 2], + output_path + 'R2s', + title, + wind=wind, + cmap='Reds' + ) + title = ( + f"R_2 (model), std={np.nanstd(ccd_maps[:, 1, 2]):.5e}\n" + + f"vmax={np.nanmax(abs(ccd_maps[:, 1, 2])):.4e}" + ) MeanShapesPlot( - ccd_maps[:, 1, 2], output_path + 'R2m', - 'R_2 (model), std=%.5e\nvmax=%.4e' % ( - np.nanstd(ccd_maps[:, 1, 2]), - np.nanmax(abs(ccd_maps[:, 1, 2]))), - wind=wind, cmap='Reds') + ccd_maps[:, 1, 2], + output_path + 'R2m', + title, + wind=wind, + cmap='Reds' + ) if auto_colorbar: wind = [0, np.nanmax(np.abs( @@ -297,44 +373,59 @@ def plot_meanshapes(starcat_path, output_path, nb_pixel, w_log, R2_res = (ccd_maps[:, 0, 2] - ccd_maps[:, 1, 2]) / ccd_maps[:, 0, 2] R2_res = R2_res[~np.isnan(R2_res)] rmse_r2 = np.sqrt(np.mean(R2_res ** 2)) - w_log.info('Bins: R2 residual RMSE: %.6f\n' % rmse_r2) + w_log.info(f"Bins: R2 residual RMSE: {rmse_r2:.6f}\n") vmax = np.nanmax( abs((ccd_maps[:, 0, 2] - ccd_maps[:, 1, 2]) / ccd_maps[:, 0, 2])) wind = [0, vmax] + std_title = np.nanstd( + (ccd_maps[:, 0, 2] - ccd_maps[:, 1, 2]) / ccd_maps[:, 0, 2] + ) + title = ( + f"āˆ†(R_2)/R_2 res, rmse={rmse_r2:.5e}\nvmax={vmax:.4e} , " + + f"std={std_title:.5e}" + ) if remove_outliers: - plot_title = "Outliers removed\nāˆ†(R_2)/R_2 res, " \ - "rmse=%.5e\nvmax=%.4e , std=%.5e" % ( - rmse_r2, vmax, - np.nanstd( - (ccd_maps[:, 0, 2] - ccd_maps[:, 1, 2]) / - ccd_maps[:, 0, 2])) - else: - plot_title = "āˆ†(R_2)/R_2 res, rmse=%.5e\nvmax=%.4e , std=%.5e" % ( - rmse_r2, vmax, - np.nanstd((ccd_maps[:, 0, 2] - ccd_maps[:, 1, 2]) / - ccd_maps[:, 0, 2])) + title = "Outliers removed\n" + title MeanShapesPlot( np.abs((ccd_maps[:, 0, 2] - ccd_maps[:, 1, 2]) / ccd_maps[:, 0, 2]), output_path + 'R2res', - plot_title, wind=wind, colorbar_ampl=colorbar_ampl, cmap='Reds') + title, + wind=wind, + colorbar_ampl=colorbar_ampl, + cmap='Reds' + ) # nstars wind = (0, np.max(ccd_maps[:, 0, 3])) + title = f"Number of stars\nTotal={np.nansum(ccd_maps[:, 0, 3]):d}" MeanShapesPlot( - ccd_maps[:, 0, 3], output_path + 'nstar', - 'Number of stars\nTotal=%d' % (np.nansum(ccd_maps[:, 0, 3])), - wind=wind, cmap='magma') + ccd_maps[:, 0, 3], + output_path + 'nstar', + title, + wind=wind, + cmap='magma' + ) # Histograms if plot_histograms: hist_bins = 50 plt.figure(figsize=(12, 6), dpi=300) - plt.hist(all_star_shapes[0, :], bins=hist_bins, range=[-0.2, 0.2], - label='stars', alpha=0.5) - plt.hist(all_psf_shapes[0, :], bins=hist_bins, range=[-0.2, 0.2], - label='PSFs', alpha=0.5) + plt.hist( + all_star_shapes[0, :], + bins=hist_bins, + range=[-0.2, 0.2], + label='stars', + alpha=0.5 + ) + plt.hist( + all_psf_shapes[0, :], + bins=hist_bins, + range=[-0.2, 0.2], + label='PSFs', + alpha=0.5 + ) plt.legend(loc='best', fontsize=16) plt.title('e1', fontsize=24) plt.savefig(output_path + 'e1_hist.png') @@ -342,19 +433,33 @@ def plot_meanshapes(starcat_path, output_path, nb_pixel, w_log, plt.figure(figsize=(12, 6), dpi=300) data_hist = all_star_shapes[0, :] - all_psf_shapes[0, :] - plt.hist(data_hist, bins=hist_bins, - range=[np.min(data_hist), np.max(data_hist)], - label='err(star - psf)', alpha=0.5) + plt.hist( + data_hist, + bins=hist_bins, + range=[np.min(data_hist), np.max(data_hist)], + label='err(star - psf)', + alpha=0.5 + ) plt.legend(loc='best', fontsize=16) plt.title('e1 err', fontsize=24) plt.savefig(output_path + 'err_e1_hist.png') plt.close() plt.figure(figsize=(12, 6), dpi=300) - plt.hist(all_star_shapes[1, :], bins=hist_bins, range=[-0.2, 0.2], - label='stars', alpha=0.5) - plt.hist(all_psf_shapes[1, :], bins=hist_bins, range=[-0.2, 0.2], - label='PSFs', alpha=0.5) + plt.hist( + all_star_shapes[1, :], + bins=hist_bins, + range=[-0.2, 0.2], + label='stars', + alpha=0.5 + ) + plt.hist( + all_psf_shapes[1, :], + bins=hist_bins, + range=[-0.2, 0.2], + label='PSFs', + alpha=0.5 + ) plt.legend(loc='best', fontsize=16) plt.title('e2', fontsize=24) plt.savefig(output_path + 'e2_hist.png') @@ -362,9 +467,13 @@ def plot_meanshapes(starcat_path, output_path, nb_pixel, w_log, plt.figure(figsize=(12, 6), dpi=300) data_hist = all_star_shapes[1, :] - all_psf_shapes[1, :] - plt.hist(data_hist, bins=hist_bins, - range=[np.min(data_hist), np.max(data_hist)], - label='err(star - psf)', alpha=0.5) + plt.hist( + data_hist, + bins=hist_bins, + range=[np.min(data_hist), np.max(data_hist)], + label='err(star - psf)', + alpha=0.5 + ) plt.legend(loc='best', fontsize=16) plt.title('e2 err', fontsize=24) plt.savefig(output_path + 'err_e2_hist.png') @@ -373,10 +482,20 @@ def plot_meanshapes(starcat_path, output_path, nb_pixel, w_log, plt.figure(figsize=(12, 6), dpi=300) mean_R2 = np.mean(all_star_shapes[2, :]) wind = [mean_R2 - 4, mean_R2 + 4] - plt.hist(all_star_shapes[2, :], bins=hist_bins, range=wind, label='stars', - alpha=0.5) - plt.hist(all_psf_shapes[2, :], bins=hist_bins, range=wind, label='PSFs', - alpha=0.5) + plt.hist( + all_star_shapes[2, :], + bins=hist_bins, + range=wind, + label='stars', + alpha=0.5 + ) + plt.hist( + all_psf_shapes[2, :], + bins=hist_bins, + range=wind, + label='PSFs', + alpha=0.5 + ) plt.legend(loc='best', fontsize=16) plt.title('R2', fontsize=24) plt.savefig(output_path + 'R2_hist.png') @@ -385,9 +504,13 @@ def plot_meanshapes(starcat_path, output_path, nb_pixel, w_log, plt.figure(figsize=(12, 6), dpi=300) data_hist = (all_star_shapes[2, :] - all_psf_shapes[2, :]) / all_star_shapes[2, :] - plt.hist(data_hist, bins=hist_bins, - range=[np.min(data_hist), np.max(data_hist)], - label='err(star - psf)/star', alpha=0.5) + plt.hist( + data_hist, + bins=hist_bins, + range=[np.min(data_hist), np.max(data_hist)], + label='err(star - psf)/star', + alpha=0.5 + ) plt.legend(loc='best', fontsize=16) plt.title('R2 err', fontsize=24) plt.savefig(output_path + 'err_R2_hist.png') @@ -399,12 +522,21 @@ def plot_meanshapes(starcat_path, output_path, nb_pixel, w_log, # Rho stats functions - -def NegDash(x_in, y_in, yerr_in, plot_name='', vertical_lines=True, - xlabel='', ylabel='', rho_nb='', - ylim=None, semilogx=False, semilogy=False, - **kwargs): - """ This function is for making plots with vertical errorbars, +def NegDash( + x_in, + y_in, + yerr_in, + plot_name='', + vertical_lines=True, + xlabel='', + ylabel='', + rho_nb='', + ylim=None, + semilogx=False, + semilogy=False, + **kwargs +): + r""" This function is for making plots with vertical errorbars, where negative values are shown in absolute value as dashed lines. The resulting plot can either be saved by specifying a file name as `plot_name', or be kept as a pyplot instance (for instance to combine @@ -487,8 +619,13 @@ class new_BaseCorrelationFunctionSysTest(BaseCorrelationFunctionSysTest): """ - def makeCatalog(self, data, config=None, use_as_k=None, - use_chip_coords=False): + def makeCatalog( + self, + data, + config=None, + use_as_k=None, + use_chip_coords=False + ): if data is None or isinstance(data, treecorr.Catalog): return data catalog_kwargs = {} @@ -543,8 +680,15 @@ class Rho1SysTest(new_BaseCorrelationFunctionSysTest): objects_list = ['star PSF'] required_quantities = [('ra', 'dec', 'g1', 'g2', 'psf_g1', 'psf_g2', 'w')] - def __call__(self, data, data2=None, random=None, random2=None, - config=None, **kwargs): + def __call__( + self, + data, + data2=None, + random=None, + random2=None, + config=None, + **kwargs + ): new_data = data.copy() new_data['g1'] = new_data['g1'] - new_data['psf_g1'] new_data['g2'] = new_data['g2'] - new_data['psf_g2'] @@ -580,8 +724,15 @@ class DESRho2SysTest(new_BaseCorrelationFunctionSysTest): objects_list = ['star PSF'] required_quantities = [('ra', 'dec', 'g1', 'g2', 'psf_g1', 'psf_g2', 'w')] - def __call__(self, data, data2=None, random=None, random2=None, - config=None, **kwargs): + def __call__( + self, + data, + data2=None, + random=None, + random2=None, + config=None, + **kwargs + ): new_data = np.rec.fromarrays([data['ra'], data['dec'], data['g1'], data['g2'], data['w']], names=['ra', 'dec', 'g1', 'g2', 'w']) @@ -625,8 +776,15 @@ class DESRho3SysTest(new_BaseCorrelationFunctionSysTest): required_quantities = [('ra', 'dec', 'sigma', 'g1', 'g2', 'psf_sigma', 'w')] - def __call__(self, data, data2=None, random=None, random2=None, - config=None, **kwargs): + def __call__( + self, + data, + data2=None, + random=None, + random2=None, + config=None, + **kwargs + ): new_data = np.rec.fromarrays( [data['ra'], data['dec'], data['g1'] * (data['sigma'] - data[ @@ -634,7 +792,8 @@ def __call__(self, data, data2=None, random=None, random2=None, data['g2'] * (data['sigma'] - data[ 'psf_sigma']) / data['sigma'], data['w']], - names=['ra', 'dec', 'g1', 'g2', 'w']) + names=['ra', 'dec', 'g1', 'g2', 'w'] + ) if data2 is not None: new_data2 = np.rec.fromarrays( [data2['ra'], data2['dec'], @@ -643,7 +802,8 @@ def __call__(self, data, data2=None, random=None, random2=None, data2['g2'] * (data2['sigma'] - data2['psf_sigma']) / data2['sigma'], data2['w']], - names=['ra', 'dec', 'g1', 'g2', 'w']) + names=['ra', 'dec', 'g1', 'g2', 'w'] + ) else: new_data2 = data2 @@ -655,7 +815,8 @@ def __call__(self, data, data2=None, random=None, random2=None, random['g2'] * (random['sigma'] - random['psf_sigma']) / random['sigma'], random['w']], - names=['ra', 'dec', 'g1', 'g2', 'w']) + names=['ra', 'dec', 'g1', 'g2', 'w'] + ) else: new_random = random @@ -667,12 +828,21 @@ def __call__(self, data, data2=None, random=None, random2=None, random2['g2'] * (random2['sigma'] - random2['psf_sigma']) / random2['sigma'], random2['w']], - names=['ra', 'dec', 'g1', 'g2', 'w']) + names=['ra', 'dec', 'g1', 'g2', 'w'] + ) else: new_random2 = random2 - return self.getCF('gg', new_data, new_data2, new_random, new_random2, - config=config, **kwargs) + + return self.getCF( + 'gg', + new_data, + new_data2, + new_random, + new_random2, + config=config, + **kwargs + ) class DESRho4SysTest(new_BaseCorrelationFunctionSysTest): @@ -686,12 +856,20 @@ class DESRho4SysTest(new_BaseCorrelationFunctionSysTest): required_quantities = [('ra', 'dec', 'g1', 'g2', 'sigma', 'psf_g1', 'psf_g2', 'psf_sigma', 'w')] - def __call__(self, data, data2=None, random=None, random2=None, - config=None, **kwargs): + def __call__( + self, + data, + data2=None, + random=None, + random2=None, + config=None, + **kwargs + ): new_data = np.rec.fromarrays( [data['ra'], data['dec'], data['g1'] - data['psf_g1'], data['g2'] - data['psf_g2'], data['w']], - names=['ra', 'dec', 'g1', 'g2', 'w']) + names=['ra', 'dec', 'g1', 'g2', 'w'] + ) if data2 is None: data2 = data new_data2 = np.rec.fromarrays( @@ -701,14 +879,16 @@ def __call__(self, data, data2=None, random=None, random2=None, data2['g2'] * (data2['sigma'] - data2['psf_sigma']) / data2[ 'sigma'], data2['w']], - names=['ra', 'dec', 'g1', 'g2', 'w']) + names=['ra', 'dec', 'g1', 'g2', 'w'] + ) if random is not None: new_random = np.rec.fromarrays( [random['ra'], random['dec'], random['g1'] - random['psf_g1'], random['g2'] - random['psf_g2'], random['w']], - names=['ra', 'dec', 'g1', 'g2', 'w']) + names=['ra', 'dec', 'g1', 'g2', 'w'] + ) else: new_random = random if random2 is None: @@ -721,11 +901,19 @@ def __call__(self, data, data2=None, random=None, random2=None, random2['g2'] * (random2['sigma'] - random2['psf_sigma']) / random2['sigma'], random2['w']], - names=['ra', 'dec', 'g1', 'g2', 'w']) + names=['ra', 'dec', 'g1', 'g2', 'w'] + ) else: new_random2 = random2 - return self.getCF('gg', new_data, new_data2, new_random, new_random2, - config=config, **kwargs) + return self.getCF( + 'gg', + new_data, + new_data2, + new_random, + new_random2, + config=config, + **kwargs + ) class DESRho5SysTest(new_BaseCorrelationFunctionSysTest): @@ -741,8 +929,15 @@ class DESRho5SysTest(new_BaseCorrelationFunctionSysTest): required_quantities = [('ra', 'dec', 'sigma', 'g1', 'g2', 'psf_sigma', 'w')] - def __call__(self, data, data2=None, random=None, random2=None, - config=None, **kwargs): + def __call__( + self, + data, + data2=None, + random=None, + random2=None, + config=None, + **kwargs + ): new_data = np.rec.fromarrays([data['ra'], data['dec'], data['g1'], data['g2'], data['w']], names=['ra', 'dec', 'g1', 'g2', 'w']) @@ -772,16 +967,28 @@ def __call__(self, data, data2=None, random=None, random2=None, random2['g2'] * (random2['sigma'] - random2['psf_sigma']) / random2['sigma'], random2['w']], - names=['ra', 'dec', 'g1', 'g2', 'w']) + names=['ra', 'dec', 'g1', 'g2', 'w'] + ) else: new_random2 = random2 - return self.getCF('gg', new_data, new_data2, new_random, new_random2, - config=config, **kwargs) - - -def rho_stats(starcat_path, output_path, rho_def='HSC', - print_fun=lambda x: print(x)): + return self.getCF( + 'gg', + new_data, + new_data2, + new_random, + new_random2, + config=config, + **kwargs + ) + + +def rho_stats( + starcat_path, + output_path, + rho_def='HSC', + print_fun=lambda x: print(x) +): """Compute and plot the five rho statistics. Syntax: @@ -815,12 +1022,14 @@ def rho_stats(starcat_path, output_path, rho_def='HSC', starcat[hdu_no].data['E2_PSF_HSM'], starcat[hdu_no].data['SIGMA_PSF_HSM'] ** 2], names=['w', 'ra', 'dec', 'g1', 'g2', 'sigma', 'psf_g1', 'psf_g2', - 'psf_sigma']) + 'psf_sigma'] + ) # TreeCorr config: TreeCorrConfig = {'ra_units': 'degrees', 'dec_units': 'degrees', 'max_sep': '3e2', 'min_sep': 5e-1, 'sep_units': 'arcmin', - 'nbins': 32} + 'nbins': 32 + } # Ininitialize all 5 rho stats if rho_def == 'HSC': diff --git a/shapepipe/modules/MCCD_package/shapepipe_auxiliary_mccd.py b/shapepipe/modules/mccd_package/shapepipe_auxiliary_mccd.py similarity index 66% rename from shapepipe/modules/MCCD_package/shapepipe_auxiliary_mccd.py rename to shapepipe/modules/mccd_package/shapepipe_auxiliary_mccd.py index 1dfc343f9..d0bd5fa52 100644 --- a/shapepipe/modules/MCCD_package/shapepipe_auxiliary_mccd.py +++ b/shapepipe/modules/mccd_package/shapepipe_auxiliary_mccd.py @@ -15,22 +15,25 @@ from astropy.io import fits import galsim from shapepipe.pipeline import file_io as sc +import pprint NOT_ENOUGH_STARS = 'Not enough stars to train the model.' -def mccd_preprocessing_pipeline(input_file_list, - output_path, - input_file_position=0, - min_n_stars=20, - separator='-', - CCD_id_filter_list=None, - outlier_std_max=100., - save_masks=True, - save_name='train_star_selection', - save_extension='.fits', - verbose=True, - print_fun=None): +def mccd_preprocessing_pipeline( + input_file_list, + output_path, + input_file_position=0, + min_n_stars=20, + separator='-', + CCD_id_filter_list=None, + outlier_std_max=100., + save_masks=True, + save_name='train_star_selection', + save_extension='.fits', + verbose=True, + print_fun=None +): r"""Preprocess input catalog. Parameters @@ -96,7 +99,9 @@ def print_fun(x): print_fun('Processing dataset..') mccd_inputs = mccd.mccd_utils.MccdInputs(separator=separator) catalog_ids = mccd_inputs.proprocess_pipeline_data( - input_file_list, element_position=input_file_position) + input_file_list, + element_position=input_file_position + ) # Loop over the catalogs for it in range(catalog_ids.shape[0]): @@ -107,8 +112,16 @@ def print_fun(x): star_list, pos_list, mask_list, ccd_list, SNR_list, RA_list, \ DEC_list, _ = mccd_inputs.outlier_rejection( - star_list, pos_list, mask_list, ccd_list, SNR_list, RA_list, - DEC_list, shape_std_max=outlier_std_max, print_fun=print_fun) + star_list, + pos_list, + mask_list, + ccd_list, + SNR_list, + RA_list, + DEC_list, + shape_std_max=outlier_std_max, + print_fun=print_fun + ) mccd_star_list = [] mccd_pos_list = [] @@ -135,14 +148,18 @@ def print_fun(x): mccd_RA_list.append(RA_list[j]) mccd_DEC_list.append(DEC_list[j]) else: - msg = '''Not enough stars in catalog_id %s - ,ccd %d. Total stars = %d''' % ( - catalog_id, ccd_list[j], n_stars) + msg = ( + f"Not enough stars in catalog_id {catalog_id} " + + f",ccd {ccd_list[j]:d}. " + + f"Total stars = {n_stars:d}." + ) print_fun(msg) except Exception: - msg = '''Warning! Problem detected in - catalog_id %s ,ccd %d''' % (catalog_id, ccd_list[j]) + msg = ( + f"Warning! Problem detected in catalog_id " + + f"{catalog_id} ,ccd {ccd_list[j]:d}" + ) print_fun(msg) if mccd_pos_list: @@ -177,30 +194,54 @@ def print_fun(x): mccd_star_nb += mccd_stars.shape[0] # Save the fits file - train_dic = {'VIGNET_LIST': mccd_stars, - 'GLOB_POSITION_IMG_LIST': mccd_poss, - 'MASK_LIST': mccd_masks, 'CCD_ID_LIST': mccd_ccds, - 'SNR_WIN_LIST': mccd_SNRs, - 'RA_LIST': mccd_RAs, 'DEC_LIST': mccd_DECs} + train_dic = { + 'VIGNET_LIST': mccd_stars, + 'GLOB_POSITION_IMG_LIST': mccd_poss, + 'MASK_LIST': mccd_masks, + 'CCD_ID_LIST': mccd_ccds, + 'SNR_WIN_LIST': mccd_SNRs, + 'RA_LIST': mccd_RAs, + 'DEC_LIST': mccd_DECs + } saving_path = output_path + save_name + separator \ + catalog_id + save_extension mccd.mccd_utils.save_to_fits(train_dic, saving_path) print_fun('Finished the training dataset processing.') - print_fun('Total stars processed = %d' % mccd_star_nb) + print_fun(f"Total stars processed = {mccd_star_nb:d}") return mccd_inputs -def mccd_fit_pipeline(trainstar_path, file_number_string, mccd_parser, - output_dir, verbose, saving_name='fitted_model'): +def mccd_fit_pipeline( + trainstar_path, + file_number_string, + mccd_parser, + output_dir, + verbose, + saving_name='fitted_model', + w_log=None +): r"""Fit the MCCD model to the observations.""" # Extract the MCCD parameters from the parser mccd_inst_kw = mccd_parser.get_instance_kw() mccd_fit_kw = mccd_parser.get_fit_kw() use_SNR_weight = mccd_parser.get_extra_kw('use_SNR_weight') + # Print the model configuration so that it is saved in log files + w_log.info('MCCD configuration parameters:') + w_log.info('[INPUTS]') + inputs_dict_str = pprint.pformat({'use_SNR_weight': use_SNR_weight}) + w_log.info(inputs_dict_str) + w_log.info('[INSTANCE]') + inst_dict_str = pprint.pformat(mccd_inst_kw) + w_log.info(inst_dict_str) + w_log.info('[FIT]') + fit_dict_str = pprint.pformat(mccd_fit_kw) + w_log.info(fit_dict_str) + w_log.info('End of MCCD configuration parameters.') + # Open fits file starcat = fits.open(trainstar_path, memmap=False) @@ -213,16 +254,23 @@ def mccd_fit_pipeline(trainstar_path, file_number_string, mccd_parser, sex_thresh=-1e5, use_SNR_weight=use_SNR_weight, verbose=verbose, - saving_name=saving_name) + saving_name=saving_name + ) starcat.close() -def mccd_validation_pipeline(teststar_path, mccd_model_path, mccd_parser, - output_dir, file_number_string, w_log, - val_saving_name): +def mccd_validation_pipeline( + teststar_path, + mccd_model_path, + mccd_parser, + output_dir, + file_number_string, + w_log, + val_saving_name +): r"""Validate the MCCD trained model against a set of observations.""" - w_log.info('Validating catalog %s..' % file_number_string) + w_log.info(f"Validating catalog {file_number_string}..") # Get MCCD parameters save_extension = '.fits' @@ -236,7 +284,8 @@ def mccd_validation_pipeline(teststar_path, mccd_model_path, mccd_parser, mccd_model_path=mccd_model_path, testcat=testcat[1], **mccd_val_kw, - sex_thresh=-1e5) + sex_thresh=-1e5 + ) testcat.close() @@ -246,15 +295,23 @@ def mccd_validation_pipeline(teststar_path, mccd_model_path, mccd_parser, # Save validation dictionary to fits file mccd.mccd_utils.save_to_fits(val_dict, val_saving_path) - w_log.info('Validation catalog < %s > saved.' % val_saving_path) + w_log.info(f"Validation catalog < {val_saving_path} > saved.") else: - w_log.info('''Fitted model corresponding to catalog %s was not - found.''' % file_number_string) - - -def mccd_interpolation_pipeline(mccd_model_path, galcat_path, pos_params, - ccd_id, saving_path, get_shapes): + w_log.info( + f"Fitted model corresponding to catalog" + + f" {file_number_string} was not found." + ) + + +def mccd_interpolation_pipeline( + mccd_model_path, + galcat_path, + pos_params, + ccd_id, + saving_path, + get_shapes +): r"""Interpolate MCCD model.""" # Import MCCD model mccd_model = mccd.mccd_quickload(mccd_model_path) @@ -273,42 +330,59 @@ def mccd_interpolation_pipeline(mccd_model_path, galcat_path, pos_params, if interp_PSFs is not None: if get_shapes: - PSF_moms = [galsim.hsm.FindAdaptiveMom( - galsim.Image(psf), strict=False) for psf in interp_PSFs] - - PSF_shapes = np.array([[moms.observed_shape.g1, - moms.observed_shape.g2, - moms.moments_sigma, - int(bool(moms.error_message))] - for moms in PSF_moms]) - - shapepipe_write_output(saving_path=saving_path, - example_fits_path=galcat_path, - get_shapes=get_shapes, - interp_PSFs=interp_PSFs, - PSF_shapes=PSF_shapes) + PSF_moms = [ + galsim.hsm.FindAdaptiveMom(galsim.Image(psf), strict=False) + for psf in interp_PSFs + ] + + PSF_shapes = np.array( + [ + [moms.observed_shape.g1, + moms.observed_shape.g2, + moms.moments_sigma, + int(bool(moms.error_message))] + for moms in PSF_moms + ] + ) + + shapepipe_write_output( + saving_path=saving_path, + example_fits_path=galcat_path, + get_shapes=get_shapes, + interp_PSFs=interp_PSFs, + PSF_shapes=PSF_shapes + ) return None else: return NOT_ENOUGH_STARS -def shapepipe_write_output(saving_path, example_fits_path, get_shapes, - interp_PSFs, PSF_shapes=None): +def shapepipe_write_output( + saving_path, + example_fits_path, + get_shapes, + interp_PSFs, + PSF_shapes=None +): r"""Save interpolated PSFs dictionary to fits file. The saved files are compatible with the previous shapepipe's standard. """ - output = sc.FITSCatalog(saving_path, - open_mode=sc.BaseCatalog.OpenMode.ReadWrite, - SEx_catalog=True) + output = sc.FITSCatalog( + saving_path, + open_mode=sc.BaseCatalog.OpenMode.ReadWrite, + SEx_catalog=True + ) if get_shapes: - data = {'VIGNET': interp_PSFs, - 'E1_PSF_HSM': PSF_shapes[:, 0], - 'E2_PSF_HSM': PSF_shapes[:, 1], - 'SIGMA_PSF_HSM': PSF_shapes[:, 2], - 'FLAG_PSF_HSM': PSF_shapes[:, 3].astype(int)} + data = { + 'VIGNET': interp_PSFs, + 'E1_PSF_HSM': PSF_shapes[:, 0], + 'E2_PSF_HSM': PSF_shapes[:, 1], + 'SIGMA_PSF_HSM': PSF_shapes[:, 2], + 'FLAG_PSF_HSM': PSF_shapes[:, 3].astype(int) + } else: data = {'VIGNET': interp_PSFs} diff --git a/shapepipe/modules/mccd_plots_runner.py b/shapepipe/modules/mccd_plots_runner.py index de072d0ac..59401982f 100644 --- a/shapepipe/modules/mccd_plots_runner.py +++ b/shapepipe/modules/mccd_plots_runner.py @@ -13,7 +13,7 @@ """ from shapepipe.modules.module_decorator import module_runner -from shapepipe.modules.MCCD_package import mccd_plot_utilities as mccd_plots +from shapepipe.modules.mccd_package import mccd_plot_utilities as mccd_plots import warnings @@ -45,24 +45,33 @@ has_mpl = False -@module_runner(input_module=['mccd_merge_starcat_runner'], version='1.0', - file_pattern=['full_starcat'], - file_ext=['.fits'], numbering_scheme='-0000000', - depends=['numpy', 'mccd', 'astropy', 'matplotlib', 'stile', - 'treecorr'], - run_method='serial') -def mccd_plots_runner(input_file_list, run_dirs, file_number_string, - config, w_log): +@module_runner( + input_module=['mccd_merge_starcat_runner'], + version='1.0', + file_pattern=['full_starcat'], + file_ext=['.fits'], + numbering_scheme='-0000000', + depends=['numpy', 'mccd', 'astropy', 'matplotlib', 'stile', 'treecorr'], + run_method='serial' +) +def mccd_plots_runner( + input_file_list, + run_dirs, + file_number_string, + config, + module_config_sec, + w_log +): # Get parameters for meanshapes plots - x_nb_bins = config.getint('MCCD_PLOTS_RUNNER', 'X_GRID') - y_nb_bins = config.getint('MCCD_PLOTS_RUNNER', 'Y_GRID') - remove_outliers = config.getboolean('MCCD_PLOTS_RUNNER', 'REMOVE_OUTLIERS') - plot_meanshapes = config.getboolean('MCCD_PLOTS_RUNNER', 'PLOT_MEANSHAPES') - plot_histograms = config.getboolean('MCCD_PLOTS_RUNNER', 'PLOT_HISTOGRAMS') + x_nb_bins = config.getint(module_config_sec, 'X_GRID') + y_nb_bins = config.getint(module_config_sec, 'Y_GRID') + remove_outliers = config.getboolean(module_config_sec, 'REMOVE_OUTLIERS') + plot_meanshapes = config.getboolean(module_config_sec, 'PLOT_MEANSHAPES') + plot_histograms = config.getboolean(module_config_sec, 'PLOT_HISTOGRAMS') # Get parameters for Rho stats plots - plot_rho_stats = config.getboolean('MCCD_PLOTS_RUNNER', 'PLOT_RHO_STATS') - rho_stat_plot_style = config.get('MCCD_PLOTS_RUNNER', 'RHO_STATS_STYLE') + plot_rho_stats = config.getboolean(module_config_sec, 'PLOT_RHO_STATS') + rho_stat_plot_style = config.get(module_config_sec, 'RHO_STATS_STYLE') nb_pixel = x_nb_bins, y_nb_bins starcat_path = input_file_list[0][0] @@ -70,37 +79,49 @@ def mccd_plots_runner(input_file_list, run_dirs, file_number_string, if plot_meanshapes or plot_histograms: if has_mpl: - mccd_plots.plot_meanshapes(starcat_path, output_path, nb_pixel, - w_log, remove_outliers, - plot_meanshapes, - plot_histograms) + mccd_plots.plot_meanshapes( + starcat_path, + output_path, + nb_pixel, + w_log, + remove_outliers, + plot_meanshapes, + plot_histograms + ) else: - msg = "[!] In order to plot the Meanshapes the package " \ - "_matplotlib_ has to be correctly imported." \ - "This was not the case, so the task is" \ - "aborted. For the next time make sure the package is" \ - "installed." + msg = ( + "[!] In order to plot the Meanshapes the package " + + "_matplotlib_ has to be correctly imported. This was not" + + " the case, so the task is aborted. For the next time make" + + " sure the package is installed." + ) warnings.warn(msg) w_log.info(msg) if plot_rho_stats: if has_stile is False or has_treecorr is False: - msg = "[!] In order to calculate the Rho stats the packages " \ - "_stile_ and _treecorr_ have to be correctly imported." \ - "This was not the case, so the rho stat calculation is" \ - "aborted. For the next time make sure both of the" \ - "packages are installed." + msg = ( + "[!] In order to calculate the Rho stats the packages " + + "_stile_ and _treecorr_ have to be correctly imported." + + " This was not the case, so the rho stat calculation is" + + "aborted. For the next time make sure both of the" + + "packages are installed." + ) warnings.warn(msg) w_log.info(msg) elif rho_stat_plot_style != 'HSC' and rho_stat_plot_style != 'DEC': - msg = "The rho stat definition should be HSC or DEC." \ - "An unknown definition was used. Rho stat calculation" \ - "aborted." + msg = ( + "The rho stat definition should be HSC or DEC. An unknown" + + " definition was used. Rho stat calculation aborted." + ) warnings.warn(msg) w_log.info(msg) else: - mccd_plots.rho_stats(starcat_path, output_path, - rho_def=rho_stat_plot_style, - print_fun=lambda x: w_log.info(x)) + mccd_plots.rho_stats( + starcat_path, + output_path, + rho_def=rho_stat_plot_style, + print_fun=lambda x: w_log.info(x) + ) return None, None diff --git a/shapepipe/modules/mccd_preprocessing_runner.py b/shapepipe/modules/mccd_preprocessing_runner.py index 86fa50cd5..025099f32 100644 --- a/shapepipe/modules/mccd_preprocessing_runner.py +++ b/shapepipe/modules/mccd_preprocessing_runner.py @@ -9,21 +9,31 @@ """ from shapepipe.modules.module_decorator import module_runner -from shapepipe.modules.MCCD_package import shapepipe_auxiliary_mccd as aux_mccd +from shapepipe.modules.mccd_package import shapepipe_auxiliary_mccd\ + as aux_mccd import mccd -@module_runner(input_module=['setools_runner'], version='1.0', - file_pattern=['star_split_ratio_80', 'star_split_ratio_20'], - file_ext=['.fits', '.fits'], - depends=['numpy', 'mccd', 'galsim', 'astropy'], - run_method='serial') -def mccd_preprocessing_runner(input_file_list, run_dirs, file_number_string, - config, w_log): +@module_runner( + input_module=['setools_runner'], + version='1.0', + file_pattern=['star_split_ratio_80', 'star_split_ratio_20'], + file_ext=['.fits', '.fits'], + depends=['numpy', 'mccd', 'galsim', 'astropy'], + run_method='serial' +) +def mccd_preprocessing_runner( + input_file_list, + run_dirs, + file_number_string, + config, + module_config_sec, + w_log +): # Recover the MCCD config file and its params - config_file_path = config.getexpanded('MCCD', 'CONFIG_PATH') - mccd_mode = config.get('MCCD', 'MODE') - verbose = config.getboolean('MCCD', 'VERBOSE') + config_file_path = config.getexpanded(module_config_sec, 'CONFIG_PATH') + mccd_mode = config.get(module_config_sec, 'MODE') + verbose = config.getboolean(module_config_sec, 'VERBOSE') # Parse MCCD config file mccd_parser = mccd.auxiliary_fun.MCCDParamsParser(config_file_path) @@ -56,28 +66,32 @@ def mccd_preprocessing_runner(input_file_list, run_dirs, file_number_string, min_n_stars_list = [1] else: - raise ValueError('''MODE should be in ["FIT", "FIT_VALIDATION", - "VALIDATION"].''') + raise ValueError( + "MODE should be in ['FIT', FIT_VALIDATION', 'VALIDATION']." + ) # Use the outfile from the pipeline and ignore the output directory from # the MCCD config file # Output paths for both newly generates datasets output_mccd_path = run_dirs['output'] + '/' - [aux_mccd.mccd_preprocessing_pipeline( - input_file_list=input_file_list, - output_path=output_mccd_path, - input_file_position=_input_pos, - min_n_stars=_min_stars, - separator=separator, - CCD_id_filter_list=None, - outlier_std_max=outlier_std_max, - save_masks=False, - save_name=_save_name, - save_extension='.fits', - verbose=verbose, - print_fun=w_log.info) - for _input_pos, _save_name, _min_stars in - zip(input_file_pos_list, save_name_list, min_n_stars_list)] + [ + aux_mccd.mccd_preprocessing_pipeline( + input_file_list=input_file_list, + output_path=output_mccd_path, + input_file_position=_input_pos, + min_n_stars=_min_stars, + separator=separator, + CCD_id_filter_list=None, + outlier_std_max=outlier_std_max, + save_masks=False, + save_name=_save_name, + save_extension='.fits', + verbose=verbose, + print_fun=w_log.info + ) + for _input_pos, _save_name, _min_stars in + zip(input_file_pos_list, save_name_list, min_n_stars_list) + ] return None, None diff --git a/shapepipe/modules/mccd_val_runner.py b/shapepipe/modules/mccd_val_runner.py index e1a8828e5..a9ab9c0fb 100644 --- a/shapepipe/modules/mccd_val_runner.py +++ b/shapepipe/modules/mccd_val_runner.py @@ -9,21 +9,33 @@ """ from shapepipe.modules.module_decorator import module_runner -from shapepipe.modules.MCCD_package import shapepipe_auxiliary_mccd as aux_mccd +from shapepipe.modules.mccd_package import shapepipe_auxiliary_mccd as\ + aux_mccd import mccd -@module_runner(input_module=['mccd_preprocessing_runner'], version='1.0', - file_pattern=['fitted_model', 'test_star_selection'], - file_ext=['.npy', '.fits'], numbering_scheme='-0000000', - depends=['numpy', 'mccd', 'galsim'], run_method='parallel') -def mccd_val_runner(input_file_list, run_dirs, file_number_string, - config, w_log): +@module_runner( + input_module=['mccd_preprocessing_runner'], + version='1.0', + file_pattern=['fitted_model', 'test_star_selection'], + file_ext=['.npy', '.fits'], + numbering_scheme='-0000000', + depends=['numpy', 'mccd', 'galsim'], + run_method='parallel' +) +def mccd_val_runner( + input_file_list, + run_dirs, + file_number_string, + config, + module_config_sec, + w_log +): # Recover the MCCD config file and its params - config_file_path = config.getexpanded('MCCD', 'CONFIG_PATH') - mccd_mode = config.get('MCCD', 'MODE') - verbose = config.getboolean('MCCD', 'VERBOSE') + config_file_path = config.getexpanded(module_config_sec, 'CONFIG_PATH') + mccd_mode = config.get(module_config_sec, 'MODE') + verbose = config.getboolean(module_config_sec, 'VERBOSE') # Parse MCCD config file mccd_parser = mccd.auxiliary_fun.MCCDParamsParser(config_file_path) @@ -47,10 +59,12 @@ def mccd_val_runner(input_file_list, run_dirs, file_number_string, output_dir=output_dir, file_number_string=file_number_string, w_log=w_log, - val_saving_name=val_saving_name) + val_saving_name=val_saving_name + ) else: - raise ValueError('''mccd_val_runner should be called when the - MODE is "VALIDATION".''') + raise ValueError( + "mccd_val_runner should be called when the MODE is 'VALIDATION'." + ) return None, None