added description in functions

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()}

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')
 
 

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):