Skip to content

Commit

Permalink
Stop the decomposition procedure once find the templates coverage is …
Browse files Browse the repository at this point in the history 
…not enough.
  • Loading branch information
@WenkeRen
WenkeRen committed Jul 11, 2024
1 parent 5e431d4 commit 4b203da
Show file tree
Hide file tree
Showing 2 changed files with 56 additions and 35 deletions.
72 changes: 43 additions & 29 deletions src/pyqsofit/HostDecomp.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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')
Expand All @@ -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
Expand Down Expand Up @@ -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)

Expand All @@ -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)
Expand Down
19 changes: 13 additions & 6 deletions src/pyqsofit/PyQSOFit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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([])
Expand All @@ -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
Expand All @@ -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) / (
Expand Down

0 comments on commit 4b203da

Please sign in to comment.