From 4b203daefcd02607366bccc3ea6a7294d33a9884 Mon Sep 17 00:00:00 2001 From: Wenke Date: Thu, 11 Jul 2024 10:05:58 +0900 Subject: [PATCH] Stop the decomposition procedure once find the templates coverage is not enough. --- src/pyqsofit/HostDecomp.py | 72 +++++++++++++++++++++++--------------- src/pyqsofit/PyQSOFit.py | 19 ++++++---- 2 files changed, 56 insertions(+), 35 deletions(-) diff --git a/src/pyqsofit/HostDecomp.py b/src/pyqsofit/HostDecomp.py index a4ca7a8..05b41b7 100644 --- a/src/pyqsofit/HostDecomp.py +++ b/src/pyqsofit/HostDecomp.py @@ -9,6 +9,8 @@ # Description：Use the PCA data to decompose the spectra with a priori """ import os, glob +import warnings + from astropy.io import fits from numpy import genfromtxt import numpy as np @@ -279,6 +281,7 @@ def __init__(self, wave, flux, err, n_gal, n_qso, path, host_type='PCA', qso_typ self.err = err self.n_gal = n_gal self.n_qso = n_qso + self.assertion = True path2prior = os.path.join(path, 'pca/prior') path2qso = os.path.join(path, 'pca/Yip_pca_templates') @@ -297,27 +300,29 @@ def __init__(self, wave, flux, err, n_gal, n_qso, path, host_type='PCA', qso_typ wave_max = np.min([np.max(wave), np.max(self.qso_tmp.wave_qso), np.max(self.gal_tmp.wave_gal)]) ind_data = np.where((wave > wave_min) & (wave < wave_max), True, False) if np.sum(ind_data)/len(ind_data) < 0.5: - raise ValueError('The templates used for decomposition can only cover less than 50% of the original data. ' + self.assertion = False + warnings.warn('The templates used for decomposition can only cover less than 50% of the original data. ' 'Please check the settings and consider close the decomposition function.') - self.wave, self.flux, self.err = wave[ind_data], flux[ind_data], err[ind_data] - - if na_mask == True: - self.wave_fit, self.flux_fit, self.err_fit = _na_mask(self.wave, self.flux, self.err) else: - self.wave_fit, self.flux_fit, self.err_fit = self.wave, self.flux, self.err + self.wave, self.flux, self.err = wave[ind_data], flux[ind_data], err[ind_data] + + if na_mask == True: + self.wave_fit, self.flux_fit, self.err_fit = _na_mask(self.wave, self.flux, self.err) + else: + self.wave_fit, self.flux_fit, self.err_fit = self.wave, self.flux, self.err - self.qso_datacube = self.qso_tmp.interp_data(self.wave_fit) - self.gal_datacube = self.gal_tmp.interp_data(self.wave_fit) + self.qso_datacube = self.qso_tmp.interp_data(self.wave_fit) + self.gal_datacube = self.gal_tmp.interp_data(self.wave_fit) - self.n_qso = self.qso_tmp.n_template - self.n_gal = self.gal_tmp.n_template + self.n_qso = self.qso_tmp.n_template + self.n_gal = self.gal_tmp.n_template - if self.n_qso != n_qso: - raise ValueError(f'The number of qso template is not correct. ' - f'You ask for {n_qso} template, while the maximum number of qso template is {self.n_qso}') - if self.n_gal != n_gal: - raise ValueError(f'The number of gal template is not correct. ' - f'You ask for {n_gal} template, while the maximum number of gal template is {self.n_gal}') + if self.n_qso != n_qso: + raise ValueError(f'The number of qso template is not correct. ' + f'You ask for {n_qso} template, while the maximum number of qso template is {self.n_qso}') + if self.n_gal != n_gal: + raise ValueError(f'The number of gal template is not correct. ' + f'You ask for {n_gal} template, while the maximum number of gal template is {self.n_gal}') def auto_decomp(self): flux_temp = np.vstack((self.qso_datacube, self.gal_datacube)).T @@ -384,6 +389,13 @@ def __init__(self, wave, flux, err, n_gal, n_qso, path, host_type='PCA', qso_typ QSO_prior_name = f'QSO_pp_prior_{qso_type}.csv' self.fh_ini_list = fh_ini_list + self.wave = wave + self.flux = flux + self.err = err + self.n_gal = n_gal + self.n_qso = n_qso + self.assertion = True + self.qso_tmp = QSO_PCA(path2qso, n_qso, template_name=qso_type) self.gal_tmp = host_template(n_template=n_gal, template_path=path2host, template_type=host_type) @@ -392,24 +404,26 @@ def __init__(self, wave, flux, err, n_gal, n_qso, path, host_type='PCA', qso_typ wave_max = np.min([np.max(wave), np.max(self.qso_tmp.wave_qso), np.max(self.gal_tmp.wave_gal)]) ind_data = np.where((wave > wave_min) & (wave < wave_max), True, False) if np.sum(ind_data)/len(ind_data) < 0.5: - raise ValueError('The templates used for decomposition can only cover less than 50% of the original data. ' + warnings.warn('The templates used for decomposition can only cover less than 50% of the original data. ' 'Please check the settings and consider close the decomposition function.') - self.wave, self.flux, self.err = wave[ind_data], flux[ind_data], err[ind_data] - - if na_mask == True: - self.wave_fit, self.flux_fit, self.err_fit = _na_mask(self.wave, self.flux, self.err) + self.assertion = False else: - self.wave_fit, self.flux_fit, self.err_fit = self.wave, self.flux, self.err + self.wave, self.flux, self.err = wave[ind_data], flux[ind_data], err[ind_data] + + if na_mask == True: + self.wave_fit, self.flux_fit, self.err_fit = _na_mask(self.wave, self.flux, self.err) + else: + self.wave_fit, self.flux_fit, self.err_fit = self.wave, self.flux, self.err - self.qso_datacube = self.qso_tmp.interp_data(self.wave_fit) - self.gal_datacube = self.gal_tmp.interp_data(self.wave_fit) - self.n_qso = self.qso_tmp.n_template - self.n_gal = self.gal_tmp.n_template + self.qso_datacube = self.qso_tmp.interp_data(self.wave_fit) + self.gal_datacube = self.gal_tmp.interp_data(self.wave_fit) + self.n_qso = self.qso_tmp.n_template + self.n_gal = self.gal_tmp.n_template - self.qso_prior = np.array([[0, 1]] * self.n_qso) - self.gal_prior = np.array([[0, 1]] * self.n_gal) + self.qso_prior = np.array([[0, 1]] * self.n_qso) + self.gal_prior = np.array([[0, 1]] * self.n_gal) - self.read_prior(path2prior, QSO_prior_name, GAL_prior_name) + self.read_prior(path2prior, QSO_prior_name, GAL_prior_name) def read_prior(self, prior_loc, QSO_prior_name='QSO_pp_prior.csv', GAL_prior_name='GAL_pp_prior.csv'): self.qso_prior = self._read_prior(os.path.join(prior_loc, QSO_prior_name), self.n_qso) diff --git a/src/pyqsofit/PyQSOFit.py b/src/pyqsofit/PyQSOFit.py index 35c523d..5d9cbcb 100644 --- a/src/pyqsofit/PyQSOFit.py +++ b/src/pyqsofit/PyQSOFit.py @@ -241,7 +241,7 @@ def Fit(self, name=None, nsmooth=1, and_mask=False, or_mask=False, reject_badpix '{1}ZBIN{2}' where 1 is the luminosity bin from one of [A, B, C, D] and 2 is the redshift bin from one of [1, 2, 3, 4, 5]. Yip et al. (2004) built a series sets of quasar PCA templates based on different redshift and absolute i-band magnitude subsamples. Check https://doi.org/10.1086/425626 for more detail. If the - host_prior is set to True, then only 'global' and 'CZBIN1' is supported. + host_prior is set to True, then only 'DZBIN1' and 'CZBIN1' is supported. npca_gal: int, optional the number of galaxy PCA components applied. It only works when decompose_host is True. The default is 5, @@ -770,7 +770,7 @@ def decompose_host_qso(self, wave, flux, err, path): """Decompose the host galaxy from QSO""" # Initialize default values self.host = np.zeros(len(wave)) - self.decomposed = False + self.decomposed = True self.host_result = np.array([]) self.host_result_type = np.array([]) self.host_result_name = np.array([]) @@ -779,12 +779,20 @@ def decompose_host_qso(self, wave, flux, err, path): prior_fitter = Prior_decomp(self.wave, self.flux, self.err, self.npca_gal, self.npca_qso, path, host_type=self.host_type, qso_type=self.qso_type, na_mask=self.decomp_na_mask) - datacube, frac_host_4200, frac_host_5100, qso_par, gal_par = prior_fitter.auto_decomp(self.host_prior_scale) + if prior_fitter.assertion is True: + datacube, frac_host_4200, frac_host_5100, qso_par, gal_par = prior_fitter.auto_decomp(self.host_prior_scale) + else: + self.decomposed = False + return self.wave, self.flux, self.err else: linear_fitter = Linear_decomp(self.wave, self.flux, self.err, self.npca_gal, self.npca_qso, path, host_type=self.host_type, qso_type=self.qso_type, na_mask=self.decomp_na_mask) - datacube, frac_host_4200, frac_host_5100, qso_par, gal_par = linear_fitter.auto_decomp() + if linear_fitter.assertion is True: + datacube, frac_host_4200, frac_host_5100, qso_par, gal_par = linear_fitter.auto_decomp() + else: + self.decomposed = False + return self.wave, self.flux, self.err # for some negative host template, we do not do the decomposition # not apply anymore # For a few cases, the host template is too weak that the host spectra (data - qso) would be mostly negative @@ -793,11 +801,10 @@ def decompose_host_qso(self, wave, flux, err, path): host_spec = datacube[1, :] - datacube[4, :] if np.sum(np.where(datacube[3, :] < 0, True, False) | np.where(datacube[4, :] < 0, True, False)) > 0.1 * \ datacube.shape[1] or np.median(datacube[3, :]) < 0.01 * flux_level or np.median(host_spec) < 0: + self.decomposed = False if self.verbose: print('Got negative host galaxy / QSO flux over 10% of coverage, decomposition is not applied!') else: - self.decomposed = True - del self.wave, self.flux, self.err self.wave = datacube[0, :] rchi2_decomp = np.sum((datacube[1, :] - datacube[4, :] - datacube[3, :]) ** 2 / datacube[2, :] ** 2) / (