diff --git a/py/redrock/archetypes.py b/py/redrock/archetypes.py index 43ce4cb1..23613459 100644 --- a/py/redrock/archetypes.py +++ b/py/redrock/archetypes.py @@ -25,8 +25,6 @@ ## Loading some global data once as it will be used many times archetype_galaxies = params_for_all_galaxies() ## Dictionary for ELG, LRG, BGS ## Dictionary of galaxy properties of archetypes -rrarchs_params = return_galaxy_archetype_properties('/global/cfs/cdirs/desi/users/abhijeet/new-archetypes/rrarchetype-galaxy.fits') - class Archetype(): """Class to store all different archetypes from the same spectype. @@ -53,6 +51,7 @@ def __init__(self, filename): self._version = hdr['VERSION'] self.wave = np.asarray(hdr['CRVAL1'] + hdr['CDELT1']*np.arange(self.flux.shape[1])) + self.filename = filename if hdr['LOGLAM']: self.wave = 10**self.wave @@ -102,6 +101,8 @@ def get_best_archetype(self,spectra,weights,flux,wflux,dwave,z,legendre, nearest dwave (dic): dictionary of wavelength grids z (float): best redshift legendre (dic): legendre polynomial + nearest_nbh (bool): nearest-neighbour approach on archetypes + n_nbh (int): number of nearest neighbours taken into account (if nearest_nbh is True) Returns: chi2 (float): chi2 of best archetype @@ -128,7 +129,8 @@ def get_best_archetype(self,spectra,weights,flux,wflux,dwave,z,legendre, nearest spectype = self._full_type[iBest].split(':::')[0] #redrock best spectype subtype = self._full_type[iBest].split(':::')[1].split('_')[0] #redrock best subtype if spectype=='GALAXY': - new_arch, gal_inds = return_N_nearest_archetypes_from_synthetic_spectra(arch_id=iBest, archetype_data=rrarchs_params, gal_data=archetype_galaxies[subtype], n_nbh=n_nbh, ret_wave=False) + rrarchs_params = return_galaxy_archetype_properties(self.filename) + new_arch, gal_inds = return_N_nearest_archetypes_from_synthetic_spectra(arch_id=iBest, archetype_data=rrarchs_params, gal_data=archetype_galaxies[subtype], n_nbh=n_nbh) new_arch = new_arch.astype('float64') # to let it work with trapz_rebin and Numba for i in range(n_nbh): binned = {hs:trapz_rebin(self.wave*(1.+z), new_arch[i], wave) for hs, wave in dwave.items()} diff --git a/py/redrock/external/desi.py b/py/redrock/external/desi.py index a4a28b0f..c23049e2 100644 --- a/py/redrock/external/desi.py +++ b/py/redrock/external/desi.py @@ -560,7 +560,7 @@ def rrdesi(options=None, comm=None): parser.add_argument("--nearest_nbh", default=False, action="store_true", required=False, help="Will apply the nearest neighbour approach on archetypes") - parser.add_argument("-n_nbh", type=int, default=9, + parser.add_argument("-n_nbh", "--n_nearest", type=int, default=9, required=False, help="if nearest_nbh True, N-nearest neighbours taken into account (default is 9)") parser.add_argument("-d", "--details", type=str, default=None, @@ -691,7 +691,7 @@ def rrdesi(options=None, comm=None): if os.path.isfile(args.archetypes): file_check = args.archetypes if os.path.isdir(args.archetypes): - file_check = args.archetypes+'rrarchetype-galaxy.fits' + file_check = args.archetypes+'rrarchetype-galaxy.fits' #must have the same naming convention try: fits.open(file_check)[2].data ## HDU2 of galaxy-archetype must contain properties except IndexError: @@ -702,9 +702,9 @@ def rrdesi(options=None, comm=None): else: sys.exit(1) print('\n=================\n') - print('Archetype with galaxy properties is provided\n') - print('Nearest neighbour approach is provided, so will apply the N-nearest neighbour approach on Redrock\n') - print('%d nearest neighbours will be used...\n'%(args.n_nbh)) + print('REDSHIFT: Archetype with galaxy properties is provided\n') + print('REDSHIFT: Nearest neighbour approach is provided, so will apply the N-nearest neighbour approach on Redrock\n') + print('REDSHIFT: %d nearest neighbours will be used...\n'%(args.n_nbh)) print('===================\n') diff --git a/py/redrock/nearest_neighbours.py b/py/redrock/nearest_neighbours.py index e72e4a01..99fc5ef3 100644 --- a/py/redrock/nearest_neighbours.py +++ b/py/redrock/nearest_neighbours.py @@ -46,7 +46,7 @@ def cartesian_dist(x_ref, Y, n_nearest): inds = np.argsort(dist) return dist, inds[0:n_nearest] -def return_N_nearest_archetypes_from_synthetic_spectra(arch_id, archetype_data, gal_data, n_nbh, ret_wave=True): +def return_N_nearest_archetypes_from_synthetic_spectra(arch_id, archetype_data, gal_data, n_nbh): """ Returns the fluxes and rest-frame wavelength (same as archetypes) of N-nearest neighbours for a given Best archetypes after estimating the distance between it and all archetypes @@ -55,26 +55,17 @@ def return_N_nearest_archetypes_from_synthetic_spectra(arch_id, archetype_data, arch_id: [int]; best archetype id found by redrock in archetype mode archetype_data: [dict]; dictionary containing physical parameters for archetypes must contain: ['LOGMSTAR [M_sol]', 'LOGSSFR [yr^-1]', 'AV_ISM [mag]'] --> 3d coordinate... - gal_data: [dict]; full synthetic galaxy data dictionary, must contain:['LOGMSTAR [M_sol]', 'LOGSSFR [yr^-1]', 'AV_ISM [mag]', 'FLUX'] .... - n_nbh: [int], number of nearest neighbours you want to return - ret_wave: [bool], returns the wavelength array if True (default True) """ + param_keys = ['LOGMSTAR [M_sol]', 'LOGSSFR [yr^-1]', 'AV_ISM [mag]'] Y = np.array([gal_data[key] for key in param_keys]).T # parameters for synthetic galaxies (3d coordinates) Xref = np.array([archetype_data[key][arch_id] for key in param_keys]) # parameter for best fit redrock archetypes dist, ii = cartesian_dist(x_ref=Xref, Y=Y, n_nearest=n_nbh) arch_flux = gal_data['FLUX'][ii] - if ret_wave: - lam0 = 1228.0701754386 #starting rest frame wavelength (same as PCA templates) - npix = 19545 # number of wavelength pixels same as archetypes (5 times lower than PCA) - step_size = 0.5 # pixel size in Angstrom (5 times lower resolution than PCA templates) - temp_wave = lam0 + np.arange(npix)*step_size # Galaxy archetype wavelengths - return arch_flux, temp_wave, ii - else: - return arch_flux, ii + return arch_flux, ii def return_galaxy_archetype_properties(archetypes):