diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index df2062a8..918448df 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -96,7 +96,7 @@ jobs: fail-fast: false matrix: os: [ubuntu-latest] - python-version: ['3.10', '3.11'] + python-version: ['3.11'] steps: - name: Checkout code @@ -124,7 +124,7 @@ jobs: fail-fast: false matrix: os: [ubuntu-latest] - python-version: ['3.10', '3.11'] + python-version: ['3.11'] steps: - name: Checkout code diff --git a/README.rst b/README.rst index ffac919c..7e6037ce 100644 --- a/README.rst +++ b/README.rst @@ -9,8 +9,8 @@ :target: https://github.com/desihub/fastspecfit/actions :alt: GitHub Actions CI Status -.. |Coveralls Status| image:: https://coveralls.io/repos/desihub/fastspecfit/badge.svg - :target: https://coveralls.io/github/desihub/fastspecfit +.. |Coveralls Status| image:: https://coveralls.io/repos/desihub/fastspecfit/badge.svg?branch=main + :target: https://coveralls.io/github/desihub/fastspecfit?branch=main :alt: Test Coverage Status .. |Documentation Status| image:: https://readthedocs.org/projects/fastspecfit/badge/?version=latest diff --git a/bin/build-fsps-templates b/bin/build-fsps-templates index 387181c2..0779adfd 100755 --- a/bin/build-fsps-templates +++ b/bin/build-fsps-templates @@ -20,6 +20,15 @@ time python $HOME/code/desihub/fastspecfit/bin/build-fsps-templates --imf salpet time python $HOME/code/desihub/fastspecfit/bin/build-fsps-templates --imf chabrier time python $HOME/code/desihub/fastspecfit/bin/build-fsps-templates --imf kroupa +Note: to use different stellar libraries (e.g., MILES vs C3K), in a clean +terminal make sure you've git-updated a local checkout of the python-fsps and +that $SPS_HOME/src/sps_vars.f90 has been edited correctly and then: + +% echo SPS_HOME + /Users/ioannis/code/python-fsps/src/fsps/libfsps +% python -m pip install . --no-cache-dir --force-reinstall +% python -c "import fsps; sp = fsps.StellarPopulation(); print(sp.libraries)" + ''' import os, time, pdb import numpy as np @@ -106,8 +115,8 @@ def smooth_continuum(wave, flux, medbin=1000, smooth_window=200, return smooth -def build_templates(models, logages, agebins, imf='chabrier', - include_nebular=True): +def build_templates(models, logages, agebins=None, imf='chabrier', + sfh_ssp=False, include_nebular=True): nsed = len(models) @@ -115,8 +124,6 @@ def build_templates(models, logages, agebins, imf='chabrier', meta['age'] = 10**models['logage'] meta['zzsun'] = models['logmet'] meta['av'] = models['dust'] * 1.086 - #meta['fagn'] = models['fagn'] - #meta['qpah'] = models['qpah'] meta['mstar'] = np.zeros(nsed, 'f4') meta['sfr'] = np.zeros(nsed, 'f4') @@ -125,19 +132,32 @@ def build_templates(models, logages, agebins, imf='chabrier', print('Instantiating the StellarPopulation object...', end='') t0 = time.time() - sp = fsps.StellarPopulation(compute_vega_mags=False, - add_dust_emission=True, - add_neb_emission=True, - nebemlineinspec=include_nebular, - imf_type=imfdict[imf], # 0=Salpeter, 1=Chabrier, 2=Kroupa - dust_type=0, - dust_index=-0.7, - #sfh=0 # SSP parameters - #zcontinuous=1, - sfh=3, # tabular SFH parameters - zcontinuous=1, - ) + if sfh_ssp: + # SSP SFH + sp = fsps.StellarPopulation(compute_vega_mags=False, + add_dust_emission=True, + add_neb_emission=True, + nebemlineinspec=include_nebular, + imf_type=imfdict[imf], # 0=Salpeter, 1=Chabrier, 2=Kroupa + dust_type=0, + dust_index=-0.7, + sfh=0, # SSP parameters + zcontinuous=1, + ) + else: + # tabular SFH + sp = fsps.StellarPopulation(compute_vega_mags=False, + add_dust_emission=True, + add_neb_emission=True, + nebemlineinspec=include_nebular, + imf_type=imfdict[imf], # 0=Salpeter, 1=Chabrier, 2=Kroupa + dust_type=0, + dust_index=-0.7, + sfh=3, # tabular SFH parameters + zcontinuous=1, + ) print('...took {:.3f} sec'.format((time.time()-t0))) + print(sp.libraries) if include_nebular: print('Creating {} model spectra with nebular emission...'.format(nsed), end='') @@ -149,30 +169,29 @@ def build_templates(models, logages, agebins, imf='chabrier', sp.params['dust1'] = model['dust'] sp.params['dust2'] = model['dust'] sp.params['logzsol'] = model['logmet'] - #sp.params['fagn'] = model['fagn'] - #sp.params['duste_qpah'] = model['qpah'] - - # SSP spectrum - deprecated - #wave, flux = sp.get_spectrum(tage=10.**(logage - 9.0), peraa=True) - - # lookback time of constant SFR - agebin_indx = np.where(model['logage'] == np.float32(logages))[0] - agebin = agebins[agebin_indx, :][0] # Gyr - fspstime = agebin - agebin[0] # Gyr - tage = agebin[1] # time of observation [Gyr] - #print(tage, model['logage']) - - dt = np.diff(agebin) * 1e9 # [yr] - sfh = np.zeros_like(fspstime) + 1.0 / dt / 2 # [Msun/yr] - # force the SFR to go to zero at the edge - fspstime = np.hstack((fspstime, fspstime[-1]+1e-4)) - sfh = np.hstack((sfh, 0.0)) - - sp.set_tabular_sfh(fspstime, sfh) - - wave, flux = sp.get_spectrum(tage=tage, peraa=True) # tage in Gyr - #print(tage, sp.sfr) + if sfh_ssp: + # SSP + wave, flux = sp.get_spectrum(tage=10.**(model['logage'] - 9.), peraa=True) + else: + # lookback time of constant SFR + agebin_indx = np.where(model['logage'] == np.float32(logages))[0] + agebin = agebins[agebin_indx, :][0] # Gyr + fspstime = agebin - agebin[0] # Gyr + tage = agebin[1] # time of observation [Gyr] + #print(tage, model['logage']) + + dt = np.diff(agebin) * 1e9 # [yr] + sfh = np.zeros_like(fspstime) + 1. / dt #/ 2 # [Msun/yr] + + # force the SFR to go to zero at the edge + fspstime = np.hstack((fspstime, fspstime[-1]+1e-8)) + sfh = np.hstack((sfh, 0.)) + + sp.set_tabular_sfh(fspstime, sfh) + #print(tage, sp.sfr) + + wave, flux = sp.get_spectrum(tage=tage, peraa=True) # tage in Gyr #plt.clf() #I = (wave > 3500) * (wave < 9000) @@ -180,12 +199,10 @@ def build_templates(models, logages, agebins, imf='chabrier', #plt.savefig('junk2.png') #pdb.set_trace() - logage = model['logage'] - - lodot = 3.828 # 10^{33} erg/s - tenpc2 = (10.0 * 3.085678)**2 * 1e3 # 10^{33} cm^2 + lodot = 3.828e33 # [erg/s] + tenpc2 = (10. * 3.085678e18)**2 # [cm^2] - flux = flux * lodot / (4.0 * np.pi * tenpc2) + flux = flux * lodot / (4. * np.pi * tenpc2) # [erg/s/cm2/A/Msun at 10pc] # Resample to constant log-lambda / velocity. In the IR (starting at ~1 # micron), take every fourth sampling, to save space. @@ -210,13 +227,16 @@ def build_templates(models, logages, agebins, imf='chabrier', meta['mstar'][imodel] = sp.stellar_mass meta['sfr'][imodel] = sp.sfr + #print(tage, sp.formed_mass, sp.stellar_mass) + if sp.stellar_mass < 0: + raise ValueError('Stellar mass is negative!') #plt.clf() #I = np.where((wave > 3500) * (wave < 5600))[0] #J = np.where((newwave > 3500) * (newwave < 3600))[0] #plt.plot(wave[I], flux[I]) #plt.plot(newwave[J], fluxes[J, imodel]) - #plt.show() + #plt.savefig('junk.png') #pdb.set_trace() print('...took {:.3f} min'.format((time.time()-t0)/60)) @@ -225,31 +245,45 @@ def build_templates(models, logages, agebins, imf='chabrier', def main(args): - version = '1.2.0' - # velocity dispersion grid vdispmin = 50. - vdispmax = 425. + vdispmax = 500. dvdisp = 25. nvdisp = int(np.ceil((vdispmax - vdispmin) / dvdisp)) + 1 vdisp = np.linspace(vdispmin, vdispmax, nvdisp) - ## SSP ages - deprecated - #nages = 15 - #minlogage = 7.0 # =10 Myr - #maxlogage = 10.146 # =14 Gyr - #logages = np.linspace(minlogage, maxlogage, nages) - #logages = np.linspace(minlogage, maxlogage, nages) - - # Choose lookback time bins. - agelims = np.array([0., 0.01, 0.03, 0.1, 0.33, 1.1, 3.6, 12., 14.]) # [Gyr] - nages = len(agelims) - 1 + if args.sfh_ssp: + # SSP ages + ages = np.array([0.05, 0.02, 0.065, 0.215, 0.715, 2.4, 7.8, 13.5]) # [Gyr] + logages = np.log10(1e9 * ages) # [yr] + nages = len(ages) + #nages = 8 + #minlogage = 7.0 # =10 Myr + #maxlogage = 10.146 # =14 Gyr + #logages = np.linspace(minlogage, maxlogage, nages) + agebins = None + else: + # Choose lookback time bins. + #agelims = np.array([0., 0.03, 0.1, 0.6, 3.6, 13.]) # [Gyr] + #agelims = np.array([0., 0.01, 0.03, 0.1, 0.33, 1.1, 3.6, 12., 14.]) # [Gyr] + + # from prospect.templates.adjust_continuity_agebins + nbins = 5 + tuniv = 13.7 + tbinmax = (tuniv * 0.85) * 1e9 + lim1, lim2 = 7.4772, 8.0 + agelims = np.array([0, lim1] + np.linspace(lim2, np.log10(tbinmax), nbins-2).tolist() + [np.log10(tuniv*1e9)]) # log10(yr) + agelims = 10.**agelims / 1e9 # [Gyr] + + agebins = np.array([agelims[:-1], agelims[1:]]).T # [Gyr] + logages = np.log10(1e9*np.sum(agebins, axis=1) / 2) # mean age [log10(yr)] in each bin + + print(': '+' '.join(['{:.4f} Gyr'.format(10.**logage/1e9) for logage in logages])) - agebins = np.array([agelims[:-1], agelims[1:]]).T # [Gyr] - logages = np.log10(1e9*np.sum(agebins, axis=1) / 2) # mean age [yr] in each bin + nages = len(logages) - #logmets = np.array([-1.0, -0.3, 0.0, 0.3]) - logmets = np.array([-1.0, 0.0, 0.3]) + logmets = np.array([0.0]) # [-1.0, -0.3, 0.0, 0.3] + #logmets = np.array([-1.0, 0.0, 0.3]) # [-1.0, -0.3, 0.0, 0.3] nmets = len(logmets) zsolar = 0.019 @@ -259,39 +293,21 @@ def main(args): maxlogdust = 0.477 dusts = np.hstack((mindust, np.logspace(minlogdust, maxlogdust, ndusts-1))) - #nfagns = 3 - #minfagn = 0.0 - #minlogfagn = np.log10(0.1) - #maxlogfagn = np.log10(2.0) - #fagns = np.hstack((minfagn, np.logspace(minlogfagn, maxlogfagn, nfagns-1))) - - #nqpahs = 3 - ##nqpahs = 2 - #minqlogpah = np.log10(0.1) - #maxqlogpah = np.log10(7.0) - #qpahs = np.logspace(minqlogpah, maxqlogpah, nqpahs) - # for testing - if False: + if args.test: dusts = np.array([0.0]) logmets = [0.0] - #fagns = [0.0] ndusts = 1 nmets = 1 - #nfagns = 1 dims = (nages, nmets, ndusts) - #dims = (nages, nmets, ndusts, nfagns) models_dtype = np.dtype( [('logmet', np.float32), ('dust', np.float32), - #('fagn', np.float32), - #('qpah', np.float32), ('logage', np.float32)]) # Let's be pedantic about the procedure so we don't mess up the indexing... - #dims = (nages, nmets, ndusts, nfagns, nqpahs) models = np.zeros(dims, dtype=models_dtype) for iage, logage in enumerate(logages): @@ -301,33 +317,16 @@ def main(args): models[iage, imet, idust]['dust'] = dust models[iage, imet, idust]['logage'] = logage - #for iage, logage in enumerate(logages): - # for imet, logmet in enumerate(logmets): - # for idust, dust in enumerate(dusts): - # for ifagn, fagn in enumerate(fagns): - # models[iage, imet, idust, ifagn]['logmet'] = logmet - # models[iage, imet, idust, ifagn]['dust'] = dust - # models[iage, imet, idust, ifagn]['logage'] = logage - # models[iage, imet, idust, ifagn]['fagn'] = fagn - - #for iage, logage in enumerate(logages): - # for imet, logmet in enumerate(logmets): - # for idust, dust in enumerate(dusts): - # for ifagn, fagn in enumerate(fagns): - # for iqpah, qpah in enumerate(qpahs): - # models[iage, imet, idust, ifagn, iqpah]['logmet'] = logmet - # models[iage, imet, idust, ifagn, iqpah]['dust'] = dust - # models[iage, imet, idust, ifagn, iqpah]['logage'] = logage - # models[iage, imet, idust, ifagn, iqpah]['fagn'] = fagn - # models[iage, imet, idust, ifagn, iqpah]['qpah'] = qpah - models = models.flatten() # Build models with and without line-emission. meta, wave, flux, linewaves, linefluxes = build_templates( - models, logages, agebins, include_nebular=True, imf=args.imf) + models, logages, agebins=agebins, include_nebular=True, + imf=args.imf, sfh_ssp=args.sfh_ssp) _, _, fluxnolines, _, _ = build_templates( - models, logages, agebins, include_nebular=False, imf=args.imf) + models, logages, agebins=agebins, include_nebular=False, + imf=args.imf, sfh_ssp=args.sfh_ssp) + lineflux = flux - fluxnolines #I = (wave > 3500) * (wave < 9000) @@ -343,7 +342,6 @@ def main(args): # Select just the line-free models trimmed to the 1200-10000 A wavelength # range. I = np.where((wave > 1200) * (wave < PIXKMS_WAVESPLIT))[0] - #J = np.where(meta['fagn'] == 0)[0] vdispwave = wave[I] #nvdispmodel = len(J) @@ -359,17 +357,6 @@ def main(args): vdispflux.append(gaussian_filter1d(fluxnolines[I, :], sigma=sigma, axis=0)) vdispflux = np.stack(vdispflux, axis=-1) # [npix,nvdispmodel,nvdisp] - #vdispflux = [] - #for sigma in vdisp / PIXKMS_BLU: - # vdispflux.append(gaussian_filter1d(fluxnolines[I, :][:, J], sigma=sigma, axis=0)) - # #vdispflux.append(gaussian_filter1d(normflux, sigma=sigma, axis=0)) - #vdispflux = np.stack(vdispflux, axis=-1) # [npix,nvdispmodel,nvdisp] - - #vdispflux = [] - #for sigma in vdisp / PIXKMS_BLU: - # vdispflux.append(gaussian_filter1d(fluxnolines, sigma=sigma, axis=0)) - #vdispflux = np.stack(vdispflux, axis=-1) # [npix,nvdispmodel,nvdisp] - #K = np.where((vdispwave > 3500) * (vdispwave < 4300))[0] #plt.clf() #plt.plot(vdispwave[K], fluxnolines[I, 6][J]) @@ -378,24 +365,24 @@ def main(args): #pdb.set_trace() # Write out. - outdir = os.path.join(os.environ.get('DESI_ROOT'), 'science', 'gqp', 'templates', 'fastspecfit', version) + outdir = os.path.join(os.environ.get('DESI_ROOT'), 'science', 'gqp', 'templates', 'fastspecfit', args.version) if not os.path.isdir(outdir): os.makedirs(outdir, exist_ok=True) - outfile = os.path.join(outdir, 'ftemplates-{}-{}.fits'.format(args.imf, version)) + outfile = os.path.join(outdir, 'ftemplates-{}-{}.fits'.format(args.imf, args.version)) hduflux1 = fits.PrimaryHDU(flux) hduflux1.header['EXTNAME'] = 'FLUX' - hduflux1.header['VERSION'] = version + hduflux1.header['VERSION'] = args.version hduflux1.header['BUNIT'] = 'erg/(s cm2 Angstrom)' hduflux2 = fits.ImageHDU(lineflux) hduflux2.header['EXTNAME'] = 'LINEFLUX' - hduflux2.header['VERSION'] = version + hduflux2.header['VERSION'] = args.version hduflux2.header['BUNIT'] = 'erg/(s cm2 Angstrom)' hduflux3 = fits.ImageHDU(vdispflux) hduflux3.header['EXTNAME'] = 'VDISPFLUX' - hduflux3.header['VERSION'] = version + hduflux3.header['VERSION'] = args.version hduflux3.header['VDISPMIN'] = (vdispmin, 'minimum velocity dispersion [km/s]') hduflux3.header['VDISPMAX'] = (vdispmax, 'maximum velocity dispersion [km/s]') hduflux3.header['VDISPRES'] = (dvdisp, 'velocity dispersion spacing [km/s]') @@ -424,7 +411,7 @@ def main(args): # emission lines hduflux4 = fits.ImageHDU(linefluxes) hduflux4.header['EXTNAME'] = 'LINEFLUXES' - hduflux4.header['VERSION'] = version + hduflux4.header['VERSION'] = args.version hduflux4.header['BUNIT'] = 'erg/(s cm2 Angstrom)' hduwave3 = fits.ImageHDU(linewaves) @@ -440,7 +427,11 @@ def main(args): if __name__ == '__main__': parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) - parser.add_argument('--imf', type=str, default='chabrier', choices=['chabrier', 'salpeter', 'kroupa']) + parser.add_argument('--version', required=True, help='Version number (e.g., 1.3.0)') + parser.add_argument('--sfh-ssp', action='store_true', help='Use a SSP star formation history (SFH); otherwise a tabular SFH.') + parser.add_argument('--imf', type=str, default='chabrier', choices=['chabrier', 'salpeter', 'kroupa'], + help='Initial mass function') + parser.add_argument('--test', action='store_true', help='Generate a test set of SPS models.') args = parser.parse_args() main(args) diff --git a/bin/desi-templates b/bin/desi-templates deleted file mode 100755 index 41fa3711..00000000 --- a/bin/desi-templates +++ /dev/null @@ -1,176 +0,0 @@ -#!/usr/bin/env python -"""Build spectroscopic templates. - -""" -import pdb # for debugging -import os, sys, time -import numpy as np -import fitsio -from astropy.table import Table - -from desiutil.log import get_logger -log = get_logger() - -templatedir = os.path.join(os.getenv('DESI_ROOT'), 'users', 'ioannis', 'desi-templates') -fastspecfitdir = os.path.join(os.getenv('DESI_ROOT'), 'spectro', 'fastspecfit') - -def parse(options=None): - """Parse input arguments. - - """ - import sys, argparse - - parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) - parser.add_argument('targetclass', type=str, choices=['lrg', 'elg', 'bgs', 'qso'], - help='Target class to analyze.') - parser.add_argument('--specprod', type=str, default='denali', choices=['denali', 'cascades', 'daily'], - help='Spectroscopic production to process.') - - parser.add_argument('--mp', type=int, default=1, help='Number of multiprocessing processes per MPI rank or node.') - parser.add_argument('--minperbin', type=int, default=3, help='Minimum number of galaxies per bin.') - - parser.add_argument('--minwave', type=float, default=500.0, help='Minimum output wavelength of stacked continuum spectra.') - parser.add_argument('--maxwave', type=float, default=6e4, help='Maximum output wavelength of stacked continuum spectra.') - parser.add_argument('--min-efftime', type=float, default=None, help='Minimum effective time to keep a tile (min).') - - parser.add_argument('-o', '--outdir', default=None, type=str, help='Full path to desired output directory.') - parser.add_argument('--overwrite-parent', action='store_true', help='Overwrite the parent sample.') - parser.add_argument('--overwrite-stacks', action='store_true', help='Overwrite existing stacked spectra output files.') - parser.add_argument('--overwrite-fastspec', action='store_true', help='Overwrite existing fastspec fitting results.') - parser.add_argument('--overwrite-templates', action='store_true', help='Overwrite existing templates.') - - parser.add_argument('--empca', action='store_true', help='Build the emPCA-compatible templates.') - parser.add_argument('--qa', action='store_true', help='Build QA output for a given target class.') - - if options is None: - args = parser.parse_args() - log.info(' '.join(sys.argv)) - else: - args = parser.parse_args(options) - log.info('desi-templates {}'.format(' '.join(options))) - - return args - -def main(args=None, comm=None): - """Wrapper for building templates. - - """ - if isinstance(args, (list, tuple, type(None))): - args = parse(args) - - from fastspecfit.templates.sample import SAMPLE_PROPERTIES as props - props = props[args.targetclass] - - log.info('Working on targetclass {}'.format(args.targetclass.upper())) - - tilefile = os.path.join(templatedir, '{}-tiles.fits'.format(args.targetclass)) - samplefile = os.path.join(templatedir, '{}-sample.fits'.format(args.targetclass)) - stackfile = os.path.join(templatedir, '{}-stacks.fits'.format(args.targetclass)) - fastspecfile = os.path.join(templatedir, '{}-fastspec.fits'.format(args.targetclass)) - templatefile = os.path.join(templatedir, '{}-templates.fits'.format(args.targetclass)) - - # Optional build QA for every (existing) output files and return. - if args.qa: - from fastspecfit.templates.qa import build_all_qa - build_all_qa(args.targetclass, templatedir, tilefile=tilefile, - samplefile=samplefile, stackfile=stackfile, - fastspecfile=fastspecfile, templatefile=templatefile, - specprod=args.specprod) - - return - - # [1] Build the parent sample. - if not os.path.isfile(samplefile) or args.overwrite_parent: - log.info('Building the parent sample.') - - from fastspecfit.templates.sample import select_tiles, read_fastspecfit, select_parent_sample - - # Read the master tile file and the fastspecfit fitting results and - # select the parent sample. - tilestable = select_tiles(args.targetclass, min_efftime=args.min_efftime, - remove_vi=False, specprod=args.specprod, - outfile=tilefile) - - allphot, allspec, allmeta = read_fastspecfit( - tilestable, targetclass=args.targetclass, - fastspecfit_dir=fastspecfitdir, - specprod=args.specprod) - phot, spec, meta = select_parent_sample( - allphot, allspec, allmeta, targetclass=args.targetclass, - specprod=args.specprod, samplefile=samplefile, - zobj_minmax=props['zminmax'], verbose=True) - - # [2] Build the stacked spectra. - if not os.path.isfile(stackfile) or args.overwrite_stacks: - log.info('Building stacked spectra in bins of properties.') - - from fastspecfit.templates.templates import spectra_in_bins - from fastspecfit.templates.sample import read_parent_sample, read_tilestable - - tilestable = read_tilestable(tilefile) - phot, spec, meta = read_parent_sample(samplefile) - - # select spectra in bins of properties, looping over all tiles - #tilestable = tilestable[tilestable['TILEID'] == 80607] - spectra_in_bins(tilestable, targetclass=args.targetclass, - fastspecfit_dir=fastspecfitdir, - minperbin=args.minperbin, minwave=args.minwave, - maxwave=args.maxwave, normwave=props['normwave'], - mp=args.mp, fastphot_in_bins=True, verbose=False, - stackfile=stackfile) - - # [3] Model the stacked spectra using fastspecfit. - if not os.path.isfile(fastspecfile) or args.overwrite_fastspec: - from fastspecfit.templates.templates import fastspecfit_stacks - - log.info('Modeling stacked spectra using fastspec.') - - qadir = os.path.join(templatedir, 'qa') - fastspecfit_stacks(stackfile, mp=args.mp, fastspecfile=fastspecfile, - qadir=qadir, qaprefix=args.targetclass) - - # [4] Generate the final templates. - if not os.path.isfile(templatefile) or args.overwrite_templates: - from fastspecfit.templates.templates import build_templates - - log.info('Building final templates.') - build_templates(fastspecfile, mp=args.mp, templatefile=templatefile, - minwave=args.minwave, maxwave=args.maxwave) - - # [5] Build the empca template catalog using all the target classes. - if args.empca: - from astropy.table import vstack - from fastspecfit.templates.templates import read_templates, write_templates - - log.info('Building empca template set.') - - empcafile = os.path.join(templatedir, 'empca-templates.fits') - - modelwave = None - modelflux, modeldata, weights = [], [], [] - - for targetclass in ('bgs', 'lrg', 'elg'): - templatefile = os.path.join(templatedir, '{}-templates.fits'.format(targetclass)) - - modelwave1, modelflux1, modeldata1 = read_templates(templatefile) - if modelwave is None: - modelwave = modelwave1 - assert(np.all(modelwave == modelwave1)) - modelflux.append(modelflux1) - modeldata.append(modeldata1) - weights1 = np.zeros_like(modelflux1) + modeldata1['NOBJ'].data[:, np.newaxis] - weights1 /= np.sum(modeldata1['NOBJ'].data) - weights.append(weights1) - - modelflux = np.vstack(modelflux) - weights = np.vstack(weights) - modeldata = vstack(modeldata) - - write_templates(empcafile, modelwave, modelflux, modeldata, - weights=weights, empca=True) - - pdb.set_trace() - -if __name__ == '__main__': - main() - diff --git a/bin/fastspecfit-deredshift b/bin/fastspecfit-deredshift deleted file mode 100755 index ef7d8ce8..00000000 --- a/bin/fastspecfit-deredshift +++ /dev/null @@ -1,330 +0,0 @@ -#!/usr/bin/env python -"""De-redshifting code. - -fastspecfit-deredshift --fastspecfile ./fastspect.fits - -fastspecfit-deredshift --fastphotfile /global/cfs/cdirs/desi/spectro/fastspecfit/denali/tiles/merged/fastphot-denali-cumulative.fits \ - --fastspecfile /global/cfs/cdirs/desi/spectro/fastspecfit/denali/tiles/merged/fastspec-denali-cumulative.fits - -""" -import pdb # for debugging -import os, sys, time -import numpy as np - -from fastspecfit.util import C_LIGHT - -#from redrock.rebin import trapz_rebin -from desispec.interpolation import resample_flux -from desiutil.log import get_logger -log = get_logger() - -def _deredshift_one(args): - """Multiprocessing wrapper.""" - return deredshift_one(*args) - -def deredshift_one(coaddwave, coaddflux, coaddivar, restwave, Iwave): - """QA on one spectrum.""" - resampflux = np.zeros_like(restwave) - resampivar = np.zeros_like(restwave) - if len(Iwave) > 0: - #resampflux[Iwave] = trapz_rebin(coaddwave, coaddflux, restwave[Iwave]) - flux, ivar = resample_flux(restwave[Iwave], coaddwave, coaddflux, ivar=coaddivar) - resampflux[Iwave] = flux - resampivar[Iwave] = ivar - return resampflux, resampivar - -def write_deredshifted(wave, flux, ivar, fastspec, fastphot, metadata, - outfile, specprod, coadd_type): - """Write out. - - """ - from astropy.io import fits - - t0 = time.time() - - hduflux = fits.PrimaryHDU(flux.astype('f4')) - hduflux.header['EXTNAME'] = 'FLUX' - - hduivar = fits.ImageHDU(ivar.astype('f4')) - hduivar.header['EXTNAME'] = 'IVAR' - - hduwave = fits.ImageHDU(wave.astype('f8')) - hduwave.header['EXTNAME'] = 'WAVE' - hduwave.header['BUNIT'] = 'Angstrom' - hduwave.header['AIRORVAC'] = ('vac', 'vacuum wavelengths') - - hduspec = fits.convenience.table_to_hdu(fastspec) - hduspec.header['EXTNAME'] = 'FASTSPEC' - - hduphot = fits.convenience.table_to_hdu(fastphot) - hduphot.header['EXTNAME'] = 'FASTPHOT' - - hdumeta = fits.convenience.table_to_hdu(metadata) - hdumeta.header['EXTNAME'] = 'METADATA' - hdumeta.header['SPECPROD'] = (specprod, 'spectroscopic production name') - hdumeta.header['COADDTYP'] = (coadd_type, 'spectral coadd fitted') - - hx = fits.HDUList([hduflux, hduivar, hduwave, hduspec, hduphot, hdumeta]) - hx.writeto(outfile, overwrite=True, checksum=True) - print('Writing {} spectra to {} took {:.2f} sec'.format( - len(fastphot), outfile, time.time()-t0)) - -def parse(options=None): - """Parse input arguments. - - """ - import sys, argparse - - parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) - parser.add_argument('--specprod', type=str, default='denali', choices=['denali', 'cascades', 'daily'], - help='Spectroscopic production to process.') - parser.add_argument('--coadd-type', type=str, default='cumulative', choices=['cumulative', 'pernight', 'perexp'], - help='Type of spectral coadds corresponding to the input zbestfiles.') - parser.add_argument('--tile', default=None, type=str, nargs='*', help='Tile(s) to process.') - parser.add_argument('--night', default=None, type=str, nargs='*', help='Night(s) to process (ignored if coadd-type is cumulative).') - - parser.add_argument('--targetids', type=str, default=None, help='Comma-separated list of target IDs to process.') - parser.add_argument('-n', '--ntargets', type=int, help='Number of targets to process in each file.') - parser.add_argument('--firsttarget', type=int, default=0, help='Index of first object to to process in each file (0-indexed).') - parser.add_argument('--mp', type=int, default=1, help='Number of multiprocessing processes per MPI rank or node.') - - parser.add_argument('-o', '--outdir', default=None, type=str, help='Full path to desired output directory.') - - parser.add_argument('--fastphotfile', default=None, type=str, help='Full path to fastphot fitting results.') - parser.add_argument('--fastspecfile', default=None, type=str, help='Full path to fastphot fitting results.') - - parser.add_argument('--overwrite', action='store_true', help='Overwrite any existing output files.') - - if options is None: - args = parser.parse_args() - log.info(' '.join(sys.argv)) - else: - args = parser.parse_args(options) - log.info('fastspecfit-html {}'.format(' '.join(options))) - - return args - -def main(args=None, comm=None): - """Deredshift the spectra by target class. - - """ - from astropy.table import Table - from fastspecfit.continuum import ContinuumFit - from fastspecfit.emlines import EMLineFit - from fastspecfit.io import DESISpectra, read_fastspecfit - - if isinstance(args, (list, tuple, type(None))): - args = parse(args) - - # Read the fitting results and get all the unique targetids. - if (args.fastphotfile is None or args.fastspecfile is None): - log.warning('Must provide both --fastphotfile or --fastspecfile.') - return - - fastspec, metadata, specprod, coadd_type = read_fastspecfit(args.fastspecfile) - fastphot, metaphot, _specprod, _coadd_type = read_fastspecfit(args.fastphotfile, fastphot=True) - if (specprod != _specprod) or (coadd_type != _coadd_type): - log.warning('Mis-matching specprod or coadd_type in fastspec vs fastphot fitting results!') - return - assert(np.all(fastspec['TARGETID'] == fastphot['TARGETID'])) - assert(np.all(fastspec['TARGETID'] == metaphot['TARGETID'])) - assert(np.all(metadata['TILEID'] == metaphot['TILEID'])) - - #from astropy.table import join - ## temporarily add TILEID to the tables so we can use it to join, - ## together with TARGETID and, optionally, NIGHT. - #joinkeys = ['TARGETID', 'TILEID'] - #if 'NIGHT' in fastspec.colnames: - # joinkeys += 'NIGHT' - #fastspec['TILEID'] = metaspec['TILEID'] - #fastphot['TILEID'] = metaphot['TILEID'] - # - #fastfit = join(fastspec, fastphot, join_type='outer', table_names=['SPEC', 'PHOT'], keys=joinkeys) - #metadata = join(metaspec, metaphot, join_type='outer', table_names=['SPEC', 'PHOT'], keys=joinkeys) - #assert(np.all(fastfit['TARGETID'] == metadata['TARGETID'])) - #assert(np.all(fastfit['TILEID'] == metadata['TILEID'])) - - # optionally trim to a particular tile and/or night - def _select_tiles_nights_targets(fastfit, metadata, tiles=None, nights=None): - if fastfit is None or metadata is None: - return fastfit, metadata - keep = np.ones(len(fastfit), bool) - if tiles: - tilekeep = np.zeros(len(fastfit), bool) - for tile in tiles: - tilekeep = np.logical_or(tilekeep, metadata['TILEID'].astype(str) == tile) - keep = np.logical_and(keep, tilekeep) - log.info('Keeping {} objects from tile(s) {}'.format(len(fastfit), ','.join(tiles))) - if nights and 'NIGHT' in metadata: - nightkeep = np.zeros(len(fastfit), bool) - for night in nights: - nightkeep = np.logical_or(nightkeep, metadata['NIGHT'].astype(str) == night) - keep = np.logical_and(keep, nightkeep) - log.info('Keeping {} objects from night(s) {}'.format(len(fastfit), ','.join(nights))) - return fastfit[keep], metadata[keep] - - fastspec, metadata = _select_tiles_nights_targets( - fastspec, metadata, tiles=args.tile, nights=args.night) - fastphot, metaphot = _select_tiles_nights_targets( - fastphot, metaphot, tiles=args.tile, nights=args.night) - - # parse the targetids optional input - if args.targetids: - targetids = [int(x) for x in args.targetids.split(',')] - - keep = np.where(np.isin(fastfit['TARGETID'], targetids))[0] - if len(keep) == 0: - log.warning('No matching targetids found!') - return - fastfit = fastfit[keep] - metadata = metadata[keep] - - if args.ntargets is not None: - keep = np.arange(args.ntargets) + args.firsttarget - log.info('Keeping {} targets.'.format(args.ntargets)) - fastspec = fastspec[keep] - fastphot = fastphot[keep] - metadata = metadata[keep] - - zbestdir = os.path.join(os.getenv('DESI_SPECTRO_REDUX'), specprod, 'tiles') - if args.outdir is None: - outdir = os.path.join(os.getenv('FASTSPECFIT_DATA'), specprod, 'tiles', 'deredshifted') - else: - outdir = args.outdir - if not os.path.isdir(outdir): - os.makedirs(outdir, exist_ok=True) - - # Initialize the continuum- and emission-line fitting classes. - t0 = time.time() - CFit = ContinuumFit() - EMFit = EMLineFit() - Spec = DESISpectra() - log.info('Initializing the classes took: {:.2f} sec'.format(time.time()-t0)) - - # read the tile info file for this production - tilefile = os.path.join(os.getenv('DESI_SPECTRO_REDUX'), specprod, 'tiles-{}.csv'.format(specprod)) - tileinfo = Table.read(tilefile) - tileinfo = tileinfo[np.isin(tileinfo['TILEID'], np.unique(metadata['TILEID']))] - tileinfo = tileinfo[np.argsort(tileinfo['TILEID'])] - log.info('Read survey info for {} tiles'.format(len(tileinfo))) - - pdb.set_trace() - - # rest-frame resampling parameters - obswave_min, obswave_max = 3600.0, 9800.0 - zmax = {'ELG': 1.6, 'LRG': 1.3, 'QSO': 3.5, 'BGS_ANY': 0.7} - - pixkms = 20.0 # pixel size [km/s]Traceback (most recent call last): - dlogwave = pixkms / C_LIGHT / np.log(10) # pixel size [log-lambda] - - # For each unique tile, select targets of each class and build the - # rest-frame spectra. - targetclasses = ['BGS_ANY', 'ELG', 'LRG', 'QSO']#, 'MWS_ANY' - - tall = time.time() - for tinfo in tileinfo: - tile = tinfo['TILEID'] - - if tinfo['SURVEY'].upper() == 'SV1': - from desitarget.sv1.sv1_targetmask import desi_mask#, bgs_mask, mws_mask - desibit = 'SV1_DESI_TARGET' - bgsbit = 'SV1_BGS_TARGET' - mwsbit = 'SV1_MWS_TARGET' - elif tinfo['SURVEY'].upper() == 'SV2': - from desitarget.sv2.sv2_targetmask import desi_mask#, bgs_mask, mws_mask - desibit = 'SV2_DESI_TARGET' - bgsbit = 'SV2_BGS_TARGET' - mwsbit = 'SV2_MWS_TARGET' - else: - NotImplementedError - - for targetclass in targetclasses: - - nicetargetclass = targetclass.replace('BGS_ANY', 'BGS').lower() - outfile = os.path.join(outdir, '{}-{}-restflux.fits'.format(nicetargetclass, tile)) - if os.path.isfile(outfile) and not args.overwrite: - log.info('Output file {} exists; skipping.'.format(outfile)) - continue - - targintile = np.where( - (tile == metadata['TILEID']) * - (metadata['Z'] > 0.01) * # minimum redshift! - #(metadata['Z_SPEC'] > 1.0) * - #(metadata['Z_SPEC'] < 1.1) * - metadata[desibit] & desi_mask.mask(targetclass) != 0)[0] - - log.info('Working on {} {}s on tile {}'.format(len(targintile), targetclass, tile)) - - if len(targintile) == 0: - log.info('No good {} targets in tile {}...moving on.'.format(targetclass, tile)) - continue - #targintile = targintile[:22] - - targ_fastspec = fastspec[targintile] - targ_fastphot = fastphot[targintile] - targ_metadata = metadata[targintile] - nobj = len(targ_metadata) - - # Construct the zbestfiles filenames based on the input (only - # supports coadd_type=='cumulative'). - allpetals = targ_metadata['FIBER'].data // 500 - thrunights = targ_metadata['THRUNIGHT'].astype(str).data - targetids = targ_metadata['TARGETID'] - assert(len(np.unique(thrunights)) == 1) - - zbestfiles = [] - for petal in set(allpetals): - indx = np.where((petal == allpetals))[0] - zbestfile = os.path.join(zbestdir, 'cumulative', str(tile), thrunights[indx[0]], 'zbest-{}-{}-thru{}.fits'.format( - petal, tile, thrunights[indx[0]])) - zbestfiles.append(zbestfile) - - Spec.select(zbestfiles=zbestfiles, specprod=specprod, - coadd_type='cumulative', targetids=targetids) - #dataphot = Spec.read_and_unpack(CFit, fastphot=True, synthphot=False) - dataspec = Spec.read_and_unpack(CFit, fastphot=False, synthphot=False, - remember_coadd=True) - assert(len(dataspec) == nobj) - - # now resample onto a constant log-lambda binning scale - restwave = 10**np.arange(np.log10(obswave_min / (1+zmax[targetclass])), np.log10(obswave_max), dlogwave) - restwave_min = np.min(restwave) - restwave_max = np.max(restwave) - npix = len(restwave) - - pad = 2 # [angstroms] - mpargs = [] - for iobj in np.arange(nobj): - oneplusz = 1 + dataspec[iobj]['zredrock'] - coaddwave = dataspec[iobj]['coadd_wave'] / oneplusz - coaddflux = dataspec[iobj]['coadd_flux'] * oneplusz - coaddivar = dataspec[iobj]['coadd_ivar'] / oneplusz**2 - Iwave = np.where((restwave > np.min(coaddwave+pad)) * (restwave < np.max(coaddwave-pad)))[0] - mpargs.append((coaddwave, coaddflux, coaddivar, restwave, Iwave)) - - t0 = time.time() - if args.mp > 1: - import multiprocessing - with multiprocessing.Pool(args.mp) as P: - res = P.map(_deredshift_one, mpargs) - else: - res = [deredshift_one(*_mpargs) for _mpargs in mpargs] - - res = list(zip(*res)) - restflux = np.vstack(res[0]).astype('f4') - restivar = np.vstack(res[1]).astype('f4') - log.info('De-redshifting done in {:.2f} sec'.format(time.time()-t0)) - - #import matplotlib.pyplot as plt - #plt.plot(restwave, np.sum(restflux, axis=0)) ; plt.xlim(3000, 4200) ; plt.savefig('test.png') - #pdb.set_trace() - - write_deredshifted(restwave, restflux, restivar, targ_fastspec, - targ_fastphot, targ_metadata, outfile, specprod, - coadd_type) - - log.info('De-redshifting everything took: {:.2f} sec'.format(time.time()-tall)) - -if __name__ == '__main__': - main() - diff --git a/bin/mergeall b/bin/mergeall deleted file mode 100755 index 639bc6d4..00000000 --- a/bin/mergeall +++ /dev/null @@ -1,23 +0,0 @@ -#! /bin/bash -# Merge individual files in preparation of a data release. - -mpi-fastspecfit --survey main --program dark --mp 32 --merge --overwrite -mpi-fastspecfit --survey main --program dark --mp 32 --merge --overwrite --fastphot -mpi-fastspecfit --survey main --program bright --mp 32 --merge --overwrite -mpi-fastspecfit --survey main --program bright --mp 32 --merge --overwrite --fastphot - -mpi-fastspecfit --survey sv1 --program dark --mp 32 --merge --overwrite -mpi-fastspecfit --survey sv1 --program dark --mp 32 --merge --overwrite --fastphot -mpi-fastspecfit --survey sv1 --program bright --mp 32 --merge --overwrite -mpi-fastspecfit --survey sv1 --program bright --mp 32 --merge --overwrite --fastphot - -mpi-fastspecfit --survey sv2 --program dark --mp 32 --merge --overwrite -mpi-fastspecfit --survey sv2 --program dark --mp 32 --merge --overwrite --fastphot -mpi-fastspecfit --survey sv2 --program bright --mp 32 --merge --overwrite -mpi-fastspecfit --survey sv2 --program bright --mp 32 --merge --overwrite --fastphot - -mpi-fastspecfit --survey sv3 --program dark --mp 32 --merge --overwrite -mpi-fastspecfit --survey sv3 --program dark --mp 32 --merge --overwrite --fastphot -mpi-fastspecfit --survey sv3 --program bright --mp 32 --merge --overwrite -mpi-fastspecfit --survey sv3 --program bright --mp 32 --merge --overwrite --fastphot - diff --git a/bin/mpi-fastspecfit b/bin/mpi-fastspecfit index 7226f695..68453131 100755 --- a/bin/mpi-fastspecfit +++ b/bin/mpi-fastspecfit @@ -24,8 +24,8 @@ import numpy as np from desiutil.log import get_logger log = get_logger() -def run_fastspecfit(args, comm=None, fastphot=False, specprod_dir=None, - makeqa=False, outdir_data='.', outdir_html='.'): +def run_fastspecfit(args, comm=None, fastphot=False, specprod_dir=None, makeqa=False, + outdir_data='.', templates=None, templateversion=None): import sys from desispec.parallel import stdouterr_redirected @@ -42,8 +42,7 @@ def run_fastspecfit(args, comm=None, fastphot=False, specprod_dir=None, comm=comm, specprod=args.specprod, specprod_dir=specprod_dir, coadd_type=args.coadd_type, survey=args.survey, program=args.program, healpix=args.healpix, tile=args.tile, night=args.night, - makeqa=args.makeqa, mp=args.mp, - fastphot=fastphot, outdir_data=outdir_data, outdir_html=outdir_html, + makeqa=args.makeqa, mp=args.mp, fastphot=fastphot, outdir_data=outdir_data, overwrite=args.overwrite) log.info('Planning took {:.2f} sec'.format(time.time() - t0)) else: @@ -85,6 +84,12 @@ def run_fastspecfit(args, comm=None, fastphot=False, specprod_dir=None, if args.ignore_quasarnet: cmdargs += ' --ignore-quasarnet' + if args.templates: + cmdargs += f' --templates {args.templates}' + + if args.templateversion: + cmdargs += f' --templateversion {args.templateversion}' + if args.ntargets is not None: cmdargs += f' --ntargets {args.ntargets}' @@ -175,6 +180,9 @@ def main(): parser.add_argument('--fastphot', action='store_true', help='Fit the broadband photometry.') + parser.add_argument('--templateversion', type=str, default=None, help='Template version number.') + parser.add_argument('--templates', type=str, default=None, help='Optional full path and filename to the templates.') + parser.add_argument('--merge', action='store_true', help='Merge all individual catalogs (for a given survey and program) into one large file.') parser.add_argument('--mergedir', type=str, help='Output directory for merged catalogs.') parser.add_argument('--merge-suffix', type=str, help='Filename suffix for merged catalogs.') @@ -193,7 +201,6 @@ def main(): parser.add_argument('--nolog', action='store_true', help='Do not write to the log file.') parser.add_argument('--dry-run', action='store_true', help='Generate but do not run commands.') - parser.add_argument('--outdir-html', default='$PSCRATCH/fastspecfit/html', type=str, help='Base output HTML directory.') parser.add_argument('--outdir-data', default='$PSCRATCH/fastspecfit/data', type=str, help='Base output data directory.') specprod_dir = None @@ -201,7 +208,6 @@ def main(): args = parser.parse_args() outdir_data = os.path.expandvars(args.outdir_data) - outdir_html = os.path.expandvars(args.outdir_html) if args.merge or args.mergeall or args.nompi: comm = None @@ -298,12 +304,12 @@ def main(): plan(comm=comm, specprod=args.specprod, specprod_dir=specprod_dir, coadd_type=args.coadd_type, survey=args.survey, program=args.program, healpix=args.healpix, tile=args.tile, night=args.night, - makeqa=args.makeqa, mp=args.mp, - fastphot=args.fastphot, outdir_data=outdir_data, - outdir_html=outdir_html, overwrite=args.overwrite) + makeqa=args.makeqa, mp=args.mp, fastphot=args.fastphot, + outdir_data=outdir_data, overwrite=args.overwrite) else: run_fastspecfit(args, comm=comm, fastphot=args.fastphot, specprod_dir=specprod_dir, - makeqa=args.makeqa, outdir_data=outdir_data, outdir_html=outdir_html) + makeqa=args.makeqa, outdir_data=outdir_data, templates=args.templates, + templateversion=args.templateversion) if __name__ == '__main__': main() diff --git a/doc/changes.rst b/doc/changes.rst index c30a5539..6dc6b276 100644 --- a/doc/changes.rst +++ b/doc/changes.rst @@ -2,11 +2,20 @@ Change Log ========== -2.4.4 (not released yet) +2.5.1 (not released yet) ------------------------ * +2.5.0 (2024-01-13) +------------------ + +* Address stellar mass bias; more constrained broad+narrow fitting; and max + velocity dispersion increased to 475 km/s (bump template version to 1.3.0) + [`PR #166`_]. + +.. _`PR #166`: https://github.com/desihub/fastspecfit/pull/166 + 2.4.3 (2023-12-03) ------------------ diff --git a/py/fastspecfit/continuum.py b/py/fastspecfit/continuum.py index 7dff5bb3..726dfd0f 100644 --- a/py/fastspecfit/continuum.py +++ b/py/fastspecfit/continuum.py @@ -97,14 +97,14 @@ def _smooth_continuum(wave, flux, ivar, redshift, camerapix=None, medbin=175, nsig = 3 # select just strong lines - zlinewaves = linetable['restwave'] * (1 + redshift) + zlinewaves = linetable['restwave'] * (1. + redshift) inrange = (zlinewaves > np.min(wave)) * (zlinewaves < np.max(wave)) if np.sum(inrange) > 0: linetable = linetable[inrange] linetable = linetable[linetable['amp'] >= 1] if len(linetable) > 0: for oneline in linetable: - zlinewave = oneline['restwave'] * (1 + redshift) + zlinewave = oneline['restwave'] * (1. + redshift) if oneline['isbroad']: if oneline['isbalmer']: sigma = maskkms_balmer @@ -119,7 +119,7 @@ def _smooth_continuum(wave, flux, ivar, redshift, camerapix=None, medbin=175, linemask[I] = True # Special: mask Ly-a (1215 A) - zlinewave = 1215.0 * (1 + redshift) + zlinewave = 1215. * (1. + redshift) if (zlinewave > np.min(wave)) * (zlinewave < np.max(wave)): sigma = maskkms_uv * zlinewave / C_LIGHT # [km/s --> Angstrom] I = (wave >= (zlinewave - nsig*sigma)) * (wave <= (zlinewave + nsig*sigma)) @@ -266,10 +266,10 @@ def _smooth_continuum(wave, flux, ivar, redshift, camerapix=None, medbin=175, # #xx.set_xlim(3800, 4300) # #xx.set_xlim(5200, 6050) # #xx.set_xlim(7000, 9000) - # xx.set_xlim(7000, 7800) + # xx.set_xlim(8250, 8400) #for xx in ax: # xx.set_ylim(-0.2, 1.5) - zlinewaves = linetable['restwave'] * (1 + redshift) + zlinewaves = linetable['restwave'] * (1. + redshift) linenames = linetable['name'] inrange = np.where((zlinewaves > np.min(wave)) * (zlinewaves < np.max(wave)))[0] if len(inrange) > 0: @@ -411,7 +411,7 @@ def _get_linesigma(zlinewaves, init_linesigma, label='Line', ax=None): ax = [None] * 3 # [OII] doublet, [OIII] 4959,5007 - zlinewaves = np.array([3728.483, 4960.295, 5008.239]) * (1 + redshift) + zlinewaves = np.array([3728.483, 4960.295, 5008.239]) * (1. + redshift) linesigma_narrow, linesigma_narrow_snr = _get_linesigma(zlinewaves, init_linesigma_narrow, label='Forbidden', #label='[OII]+[OIII]', ax=ax[0]) @@ -429,7 +429,7 @@ def _get_linesigma(zlinewaves, init_linesigma, label='Line', ax=None): linesigma_narrow = init_linesigma_narrow # Hbeta, Halpha - zlinewaves = np.array([4862.683, 6564.613]) * (1 + redshift) + zlinewaves = np.array([4862.683, 6564.613]) * (1. + redshift) linesigma_balmer, linesigma_balmer_snr = _get_linesigma(zlinewaves, init_linesigma_balmer, label='Balmer', #label=r'H$\alpha$+H$\beta$', ax=ax[1]) @@ -449,8 +449,8 @@ def _get_linesigma(zlinewaves, init_linesigma, label='Line', ax=None): #linesigma_balmer = init_linesigma_narrow # Lya, SiIV doublet, CIV doublet, CIII], MgII doublet - zlinewaves = np.array([1215.670, 1549.4795, 2799.942]) * (1 + redshift) - #zlinewaves = np.array([1215.670, 1398.2625, 1549.4795, 1908.734, 2799.942]) * (1 + redshift) + zlinewaves = np.array([1215.670, 1549.4795, 2799.942]) * (1. + redshift) + #zlinewaves = np.array([1215.670, 1398.2625, 1549.4795, 1908.734, 2799.942]) * (1. + redshift) linesigma_uv, linesigma_uv_snr = _get_linesigma(zlinewaves, init_linesigma_uv, label='UV/Broad', ax=ax[2]) @@ -581,7 +581,7 @@ def _kcorr_and_absmag(filters_out, band_shift): nout = len(filters_out) # note the factor of 1+band_shift - lambda_out = filters_out.effective_wavelengths.value / (1 + band_shift) + lambda_out = filters_out.effective_wavelengths.value / (1. + band_shift) # Multiply by (1+z) to convert the best-fitting model to the "rest # frame" and then divide by 1+band_shift to shift it and the wavelength @@ -954,8 +954,8 @@ def get_dn4000(wave, flam, flam_ivar=None, redshift=None, rest=True, log=None): dn4000, dn4000_ivar = 0.0, 0.0 if rest is False or redshift is not None: - restwave = wave / (1 + redshift) # [Angstrom] - flam2fnu = (1 + redshift) * restwave**2 / (C_LIGHT * 1e5) # [erg/s/cm2/A-->erg/s/cm2/Hz, rest] + restwave = wave / (1. + redshift) # [Angstrom] + flam2fnu = (1. + redshift) * restwave**2 / (C_LIGHT * 1e5) # [erg/s/cm2/A-->erg/s/cm2/Hz, rest] else: restwave = wave flam2fnu = restwave**2 / (C_LIGHT * 1e5) # [erg/s/cm2/A-->erg/s/cm2/Hz, rest] @@ -1256,7 +1256,7 @@ def __init__(self, ignore_photometry=False, fphoto=None, emlinesfile=None): super(ContinuumTools, self).__init__(ignore_photometry=ignore_photometry, fphoto=fphoto) from fastspecfit.emlines import read_emlines - + self.massnorm = 1e10 # stellar mass normalization factor [Msun] self.linetable = read_emlines(emlinesfile=emlinesfile) @@ -1320,7 +1320,7 @@ def build_linemask(wave, flux, ivar, redshift=0.0, nsig=7.0, linetable=None, # Initially, mask aggressively, especially the Balmer lines. png = None #png = 'smooth.png' - #png = '/global/homes/i/ioannis/desi-users/ioannis/tmp/smooth.png' + #png = '/global/cfs/cdirs/desi/users/ioannis/tmp/smooth.png' smooth, smoothsigma = _smooth_continuum(wave, flux, ivar, redshift, maskkms_uv=5000.0, maskkms_balmer=5000.0, maskkms_narrow=500.0, linetable=linetable, emlinesfile=emlinesfile, @@ -1329,7 +1329,7 @@ def build_linemask(wave, flux, ivar, redshift=0.0, nsig=7.0, linetable=None, # Get a better estimate of the Balmer, forbidden, and UV/QSO line-widths. png = None #png = 'linesigma.png' - #png = '/global/homes/i/ioannis/desi-users/ioannis/tmp/linesigma.png' + #png = '/global/cfs/cdirs/desi/users/ioannis/tmp/linesigma.png' linesigma_narrow, linesigma_balmer, linesigma_uv, linesigma_narrow_snr, linesigma_balmer_snr, linesigma_uv_snr = \ _estimate_linesigmas(wave, flux-smooth, ivar, redshift, png=png, log=log) @@ -1338,14 +1338,14 @@ def build_linemask(wave, flux, ivar, redshift=0.0, nsig=7.0, linetable=None, linemask_strong = np.zeros_like(linemask) # True = affected by strong emission lines. linenames = linetable['name'] - zlinewaves = linetable['restwave'] * (1 + redshift) + zlinewaves = linetable['restwave'] * (1. + redshift) lineamps = linetable['amp'] isbroads = linetable['isbroad'] * (linetable['isbalmer'] == False) isbalmers = linetable['isbalmer'] * (linetable['isbroad'] == False) png = None #png = 'linemask.png' - #png = '/global/homes/i/ioannis/desi-users/ioannis/tmp/linemask.png' + #png = '/global/cfs/cdirs/desi/users/ioannis/tmp/linemask.png' snr_strong = 1.5#3.0 inrange = (zlinewaves > np.min(wave)) * (zlinewaves < np.max(wave)) @@ -1610,7 +1610,7 @@ def templates2data(self, _templateflux, _templatewave, redshift=0.0, dluminosity if redshift > 0: ztemplatewave = templatewave * (1. + redshift) T = self.full_IGM(redshift, ztemplatewave) - T *= FLUXNORM * self.massnorm * (10.0 / (1e6 * dluminosity))**2 / (1.0 + redshift) + T *= FLUXNORM * self.massnorm * (10. / (1e6 * dluminosity))**2 / (1. + redshift) ztemplateflux = templateflux * T[:, np.newaxis] else: errmsg = 'Input redshift not defined, zero, or negative!' @@ -1809,34 +1809,6 @@ def call_nnls(modelflux, flux, ivar, xparam=None, debug=False, return chi2min, xbest, xivar, bestcoeff - @staticmethod - def get_mean_property(templateinfo, physical_property, coeff, agekeep, - normalization=None, log10=False, log=None): - """Compute the mean physical properties, given a set of coefficients. - - """ - if log is None: - from desiutil.log import get_logger - log = get_logger() - - values = templateinfo[physical_property][agekeep] # account for age of the universe trimming - - if np.count_nonzero(coeff > 0) == 0: - log.warning('Coefficients are all zero!') - meanvalue = 0.0 - #raise ValueError - else: - meanvalue = values.dot(coeff) - # the coefficients include the stellar mass normalization - if physical_property != 'mstar' and physical_property != 'sfr': - meanvalue /= np.sum(coeff) - if normalization: - meanvalue /= normalization - if log10 and meanvalue > 0: - meanvalue = np.log10(meanvalue) - - return meanvalue - def continuum_fluxes(self, data, templatewave, continuum, log=None): """Compute rest-frame luminosities and observed-frame continuum fluxes. @@ -1863,7 +1835,7 @@ def continuum_fluxes(self, data, templatewave, continuum, log=None): # luminosity-based SFRs) and at the positions of strong nebular emission # lines [OII], Hbeta, [OIII], and Halpha - dfactor = (1 + redshift) * 4.0 * np.pi * (3.08567758e24 * dlum)**2 / FLUXNORM + dfactor = (1. + redshift) * 4. * np.pi * (3.08567758e24 * dlum)**2 / FLUXNORM lums = {} cwaves = [1500.0, 2800.0, 1450., 1700., 3000., 5100.] @@ -1933,7 +1905,7 @@ def kcorr_and_absmag(self, data, templatewave, continuum, snrmin=2.0, log=None): # distance modulus, luminosity distance, and redshifted wavelength array dmod = data['dmodulus'] - ztemplatewave = templatewave * (1 + redshift) + ztemplatewave = templatewave * (1. + redshift) filters_in = self.filters[data['photsys']] @@ -1948,6 +1920,7 @@ def kcorr_and_absmag(self, data, templatewave, continuum, snrmin=2.0, log=None): return kcorr, absmag, ivarabsmag, synth_absmag, synth_maggies_in + def continuum_specfit(data, result, templatecache, fphoto=None, emlinesfile=None, constrain_age=False, no_smooth_continuum=False, ignore_photometry=False, fastphot=False, log=None, debug_plots=False, verbose=False): @@ -2014,7 +1987,7 @@ def _younger_than_universe(age, tuniv, agepad=0.5): agekeep = np.arange(nsed) nage = len(agekeep) - ztemplatewave = templatecache['templatewave'] * (1 + redshift) + ztemplatewave = templatecache['templatewave'] * (1. + redshift) # Photometry-only fitting. vdisp_nominal = templatecache['vdisp_nominal'] @@ -2204,9 +2177,9 @@ def _younger_than_universe(age, tuniv, agepad=0.5): data['apercorr'] = apercorr # needed for the line-fitting - # Performing the final fit using the line-free templates in the - # spectrum (since we mask those pixels) but the photometry - # synthesized from the templates with lines. + # Perform the final fit using the line-free templates in the spectrum + # (since we mask those pixels) but the photometry synthesized from the + # templates with lines. desitemplates_nolines, _ = CTools.templates2data( input_templateflux_nolines, templatecache['templatewave'], redshift=redshift, dluminosity=data['dluminosity'], @@ -2261,7 +2234,7 @@ def _younger_than_universe(age, tuniv, agepad=0.5): dn4000, dn4000_model)) png = None - #png = 'desi-users/ioannis/tmp/junk.png' + #png = '/global/cfs/cdirs/desi/users/ioannis/tmp/junk.png' linemask = np.hstack(data['linemask']) if np.all(coeff == 0): log.warning('Continuum coefficients are all zero.') @@ -2306,27 +2279,44 @@ def _younger_than_universe(age, tuniv, agepad=0.5): lums, cfluxes = {}, {} AV, age, zzsun, logmstar, sfr = 0.0, 0.0, 0.0, 0.0, 0.0 - #AV, age, zzsun, fagn, logmstar, sfr = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 else: kcorr, absmag, ivarabsmag, _, synth_bestmaggies = CTools.kcorr_and_absmag( data, templatecache['templatewave'], sedmodel, log=log) lums, cfluxes = CTools.continuum_fluxes(data, templatecache['templatewave'], sedmodel, log=log) - AV = CTools.get_mean_property(templatecache['templateinfo'], 'av', coeff, agekeep, log=log) # [mag] - age = CTools.get_mean_property(templatecache['templateinfo'], 'age', coeff, agekeep, normalization=1e9, log=log) # [Gyr] - zzsun = CTools.get_mean_property(templatecache['templateinfo'], 'zzsun', coeff, agekeep, log10=False, log=log) # [log Zsun] - #fagn = CTools.get_mean_property(templatecache['templateinfo'], 'fagn', coeff, agekeep, log=log) - logmstar = CTools.get_mean_property(templatecache['templateinfo'], 'mstar', coeff, agekeep, - normalization=1.0/CTools.massnorm, log10=True, log=log) # [Msun] - sfr = CTools.get_mean_property(templatecache['templateinfo'], 'sfr', coeff, agekeep, - normalization=1.0/CTools.massnorm, log10=False, log=log) # [Msun/yr] + # get the SPS properties + mstars = templatecache['templateinfo']['mstar'][agekeep] # [current mass in stars, Msun] + masstot = coeff.dot(mstars) + coefftot = np.sum(coeff) + logmstar = np.log10(CTools.massnorm * masstot) + zzsun = np.log10(coeff.dot(mstars * 10.**templatecache['templateinfo']['zzsun'][agekeep]) / masstot) # mass-weighted + AV = coeff.dot(templatecache['templateinfo']['av'][agekeep]) / coefftot # luminosity-weighted [mag] + age = coeff.dot(templatecache['templateinfo']['age'][agekeep]) / coefftot / 1e9 # luminosity-weighted [Gyr] + #age = coeff.dot(mstars * templatecache['templateinfo']['age'][agekeep]) / masstot / 1e9 # mass-weighted [Gyr] + sfr = CTools.massnorm * coeff.dot(templatecache['templateinfo']['sfr'][agekeep]) # [Msun/yr] rindx = np.argmin(np.abs(CTools.absmag_filters.effective_wavelengths.value / (1.+CTools.band_shift) - 5600)) - log.info('Mstar={:.4g} Msun, {}={:.2f} mag, A(V)={:.3f}, Age={:.3f} Gyr, SFR={:.3f} Msun/yr, Z/Zsun={:.3f}'.format( - logmstar, 'M{}'.format(CTools.absmag_bands[rindx]), absmag[rindx], AV, age, sfr, zzsun)) - #log.info('Mstar={:.4g} Msun, Mr={:.2f} mag, A(V)={:.3f}, Age={:.3f} Gyr, SFR={:.3f} Msun/yr, Z/Zsun={:.3f}, fagn={:.3f}'.format( - # logmstar, absmag[np.isin(CTools.absmag_bands, 'sdss_r')][0], AV, age, sfr, zzsun, fagn)) - + log.info(f'log(M/Msun)={logmstar:.2f}, M{CTools.absmag_bands[rindx]}={absmag[rindx]:.2f} mag, A(V)={AV:.3f}, Age={age:.3f} Gyr, SFR={sfr:.3f} Msun/yr, Z/Zsun={zzsun:.3f}') + + #import matplotlib.pyplot as plt + #W = np.where(coeff > 0)[0] + #info = templatecache['templateinfo'] + #sedwave = templatecache['templatewave'] * (1. + redshift) + #I = np.where((sedwave > 3000) * (sedwave < 40000))[0] + #plt.plot(sedwave[I], sedmodel[I]) + #for ww in W: + # plt.plot(sedwave[I], coeff[ww] * sedtemplates[I, ww], + # label=f'{info["age"][ww]/1e9:.2f} Gyr, A(V)={info["av"][ww]:.2f} mag') + #plt.legend(fontsize=10) + #plt.scatter(sedphot['lambda_eff'][:4], objflam[:4], color='k', s=100, zorder=10) + #plt.scatter(sedphot['lambda_eff'][:4], sedflam.dot(coeff)[:4], color='gray', s=100, marker='x', zorder=11) + #plt.xscale('log') + #plt.yscale('log') + #plt.savefig('/global/cfs/cdirs/desi/users/ioannis/tmp/foo.png') + # + #print(info[W]) + #print(coeff[W]) + # Pack it in and return. result['COEFF'][agekeep] = coeff result['RCHI2_PHOT'] = rchi2_phot @@ -2336,7 +2326,6 @@ def _younger_than_universe(age, tuniv, agepad=0.5): result['ZZSUN'] = zzsun result['LOGMSTAR'] = logmstar result['SFR'] = sfr - #result['FAGN'] = fagn result['DN4000_MODEL'] = dn4000_model for iband, (band, shift) in enumerate(zip(CTools.absmag_bands, CTools.band_shift)): diff --git a/py/fastspecfit/emlines.py b/py/fastspecfit/emlines.py index e9e03df0..0b5fb4ab 100644 --- a/py/fastspecfit/emlines.py +++ b/py/fastspecfit/emlines.py @@ -232,27 +232,53 @@ def __init__(self, fphoto=None, emlinesfile=None, uniqueid=None, self.minsigma_balmer_broad = 250. # minimum broad-line sigma [km/s] self.minsnr_balmer_broad = minsnr_balmer_broad # minimum broad-line S/N - def build_linemodels(self, redshift, wavelims=[3000, 10000], verbose=False, strict_finalmodel=True): - """Build all the multi-parameter emission-line models we will use. - - """ - def _print_linemodel(linemodel): + def summarize_linemodel(self, linemodel): + """Simple function to summarize an input linemodel.""" + def _print(linenames): for linename in linenames: for param in ['amp', 'sigma', 'vshift']: - I = np.where(param_names == linename+'_'+param)[0] + I = np.where(self.param_names == linename+'_'+param)[0] if len(I) == 1: I = I[0] if linemodel['tiedtoparam'][I] == -1: if linemodel['fixed'][I]: - print('{:25s} is FIXED'.format(linename+'_'+param)) + print('{:25s} is NOT FITTED'.format(linename+'_'+param)) + else: + print('{:25s} untied'.format(linename+'_'+param)) else: if linemodel['fixed'][I]: print('{:25s} tied to {:25s} with factor {:.4f} and FIXED'.format( - linename+'_'+param, param_names[linemodel['tiedtoparam'][I]], linemodel['tiedfactor'][I])) + linename+'_'+param, self.param_names[linemodel['tiedtoparam'][I]], linemodel['tiedfactor'][I])) else: print('{:25s} tied to {:25s} with factor {:.4f}'.format( - linename+'_'+param, param_names[linemodel['tiedtoparam'][I]], linemodel['tiedfactor'][I])) - + linename+'_'+param, self.param_names[linemodel['tiedtoparam'][I]], linemodel['tiedfactor'][I])) + + linenames = self.fit_linetable['name'].data + + print('---------------------') + print('UV/QSO (broad) lines:') + print('---------------------') + _print(linenames[(self.fit_linetable['isbroad'] == True) * (self.fit_linetable['isbalmer'] == False)]) + print() + print('--------------------------') + print('Broad Balmer+helium lines:') + print('--------------------------') + _print(linenames[(self.fit_linetable['isbroad'] == True) * (self.fit_linetable['isbalmer'] == True)]) + print() + print('---------------------------') + print('Narrow Balmer+helium lines:') + print('---------------------------') + _print(linenames[(self.fit_linetable['isbroad'] == False) * (self.fit_linetable['isbalmer'] == True)]) + print() + print('----------------') + print('Forbidden lines:') + print('----------------') + _print(linenames[(self.fit_linetable['isbroad'] == False) * (self.fit_linetable['isbalmer'] == False)]) + + def build_linemodels(self, redshift, wavelims=[3000, 10000], verbose=False, strict_broadmodel=True): + """Build all the multi-parameter emission-line models we will use. + + """ def _fix_parameters(linemodel, verbose=False): """Set the "fixed" attribute for all the parameters in a given linemodel.""" # First loop through all tied parameters and set fixed to the @@ -315,6 +341,7 @@ def _fix_parameters(linemodel, verbose=False): initvshift = 1.0 vmaxshift_narrow = 500.0 vmaxshift_broad = 2500.0 # 3000.0 + vmaxshift_balmer_broad = 2500. minsigma_narrow = 1.0 maxsigma_narrow = 750.0 # 500.0 @@ -411,7 +438,7 @@ def _fix_parameters(linemodel, verbose=False): for param, bounds, default in zip(['amp', 'sigma', 'vshift'], [[minamp_balmer_broad, maxamp_balmer_broad], [minsigma_balmer_broad, maxsigma_balmer_broad], - [-vmaxshift_broad, +vmaxshift_broad]], + [-vmaxshift_balmer_broad, +vmaxshift_balmer_broad]], [initamp, initsigma_broad, initvshift]): linemodel_broad['initial'][param_names == linename+'_'+param] = default linemodel_broad['bounds'][param_names == linename+'_'+param] = bounds @@ -520,17 +547,28 @@ def _fix_parameters(linemodel, verbose=False): linemodel_broad['tiedfactor'][param_names == linename+'_'+param] = 1.0 linemodel_broad['tiedtoparam'][param_names == linename+'_'+param] = np.where(param_names == 'ciii_1908_'+param)[0] - # Stephanie Juneau argues that the narrow *forbidden* line-widths - # should never be fully untied, so optionally keep them tied to - # [OIII] 5007 here; vshift can be free, however, to allow for - # wavelength calibration uncertainties. - if strict_finalmodel: - # Tie all forbidden lines to [OIII] 5007; the narrow Balmer and - # helium lines are separately tied together. - if fit_linetable['isbroad'][iline] == False and fit_linetable['isbalmer'][iline] == False and linename != 'oiii_5007': - for param in ['sigma']: + # Tie all the forbidden and narrow Balmer+helium lines *except + # [OIII] 4959,5007* to [NII] 6584 when we have broad lines. The + # [OIII] doublet frequently has an outflow component, so fit it + # separately. See the discussion at + # https://github.com/desihub/fastspecfit/issues/160 + if strict_broadmodel: + if fit_linetable['isbroad'][iline] == False and linename != 'nii_6584' and linename != 'oiii_4959' and linename != 'oiii_5007': + for param in ['sigma', 'vshift']: linemodel_broad['tiedfactor'][param_names == linename+'_'+param] = 1.0 - linemodel_broad['tiedtoparam'][param_names == linename+'_'+param] = np.where(param_names == 'oiii_5007_'+param)[0] + linemodel_broad['tiedtoparam'][param_names == linename+'_'+param] = np.where(param_names == 'nii_6584_'+param)[0] + + #if fit_linetable['isbroad'][iline] == False and linename != 'oiii_5007': + # for param in ['sigma', 'vshift']: + # linemodel_broad['tiedfactor'][param_names == linename+'_'+param] = 1.0 + # linemodel_broad['tiedtoparam'][param_names == linename+'_'+param] = np.where(param_names == 'oiii_5007_'+param)[0] + + ## Tie all forbidden lines to [OIII] 5007; the narrow Balmer and + ## helium lines are separately tied together. + #if fit_linetable['isbroad'][iline] == False and fit_linetable['isbalmer'][iline] == False and linename != 'oiii_5007': + # for param in ['sigma']: + # linemodel_broad['tiedfactor'][param_names == linename+'_'+param] = 1.0 + # linemodel_broad['tiedtoparam'][param_names == linename+'_'+param] = np.where(param_names == 'oiii_5007_'+param)[0] # Finally set the initial values and bounds on the doublet ratio parameters. for param, bounds, default in zip(['mgii_doublet_ratio', 'oii_doublet_ratio', 'sii_doublet_ratio'], @@ -564,11 +602,29 @@ def _fix_parameters(linemodel, verbose=False): linemodel_nobroad['tiedfactor'][param_names == linename+'_'+param] = 1.0 linemodel_nobroad['tiedtoparam'][param_names == linename+'_'+param] = np.where(param_names == 'halpha_broad_'+param)[0] + if strict_broadmodel: + # Tie the forbidden lines to [OIII] 5007. + if fit_linetable['isbalmer'][iline] == False and fit_linetable['isbroad'][iline] == False and linename != 'oiii_5007': + for param in ['sigma', 'vshift']: + linemodel_nobroad['tiedfactor'][param_names == linename+'_'+param] = 1.0 + linemodel_nobroad['tiedtoparam'][param_names == linename+'_'+param] = np.where(param_names == 'oiii_5007_'+param)[0] + + # Tie narrow Balmer and helium lines together. + if fit_linetable['isbalmer'][iline] and fit_linetable['isbroad'][iline] == False: + if linename == 'halpha': + for param in ['sigma', 'vshift']: + linemodel_nobroad['tiedfactor'][param_names == linename+'_'+param] = 0.0 + linemodel_nobroad['tiedtoparam'][param_names == linename+'_'+param] = -1 + else: + for param in ['sigma', 'vshift']: + linemodel_nobroad['tiedfactor'][param_names == linename+'_'+param] = 1.0 + linemodel_nobroad['tiedtoparam'][param_names == linename+'_'+param] = np.where(param_names == 'halpha_'+param)[0] + #linemodel_nobroad[np.logical_and(linemodel_nobroad['fixed'] == False, linemodel_nobroad['tiedtoparam'] == -1)] _fix_parameters(linemodel_nobroad, verbose=False) assert(np.all(linemodel_nobroad['tiedtoparam'][linemodel_nobroad['tiedfactor'] != 0] != -1)) - + return linemodel_broad, linemodel_nobroad def initial_guesses_and_bounds(self, data, emlinewave, emlineflux, log): @@ -700,8 +756,8 @@ def _linemodel_to_parameters(linemodel, fit_linetable): @staticmethod def populate_linemodel(linemodel, initial_guesses, param_bounds, log): - """Population an input linemodel with initial guesses and parameter bounds, - taking into account fixed parameters. + """Populate an input linemodel with initial guesses and parameter bounds, taking + into account fixed parameters. """ # Set initial values and bounds. @@ -710,7 +766,15 @@ def populate_linemodel(linemodel, initial_guesses, param_bounds, log): if linemodel['fixed'][iparam]: linemodel['initial'][iparam] = 0.0 # always set fixed parameter to zero else: - linemodel['initial'][iparam] = initial_guesses[param] + # Make sure the initial guess for the narrow Balmer+helium + # line is smaller than the guess for the broad-line model. + if linemodel[iparam]['isbalmer'] and 'sigma' in param: + if linemodel[iparam]['isbroad']: + linemodel['initial'][iparam] = 1.1 * initial_guesses[param] + else: + linemodel['initial'][iparam] = initial_guesses[param] + else: + linemodel['initial'][iparam] = initial_guesses[param] if param in param_bounds.keys(): linemodel['bounds'][iparam] = param_bounds[param] # set the lower boundary on broad lines to be XX times the local noise @@ -921,9 +985,6 @@ def optimize(self, linemodel, emlinewave, emlineflux, weights, redshift, out_linemodel.meta['nfev'] = fit_info['nfev'] out_linemodel.meta['status'] = fit_info['status'] - #if debug: - # pdb.set_trace() - # Get the final line-amplitudes, after resampling and convolution (see # https://github.com/desihub/fastspecfit/issues/139). Some repeated code # from build_emline_model... @@ -2422,7 +2483,9 @@ def emline_specfit(data, templatecache, result, continuummodel, smooth_continuum # Build all the emission-line models for this object. linemodel_broad, linemodel_nobroad = EMFit.build_linemodels( redshift, wavelims=(np.min(emlinewave)+5, np.max(emlinewave)-5), - verbose=False, strict_finalmodel=True) + verbose=False, strict_broadmodel=True) + #EMFit.summarize_linemodel(linemodel_broad) + #EMFit.summarize_linemodel(linemodel_nobroad) # Get initial guesses on the parameters and populate the two "initial" # linemodels; the "final" linemodels will be initialized with the @@ -2430,8 +2493,8 @@ def emline_specfit(data, templatecache, result, continuummodel, smooth_continuum initial_guesses, param_bounds = EMFit.initial_guesses_and_bounds( data, emlinewave, emlineflux, log) - for linemodel in [linemodel_broad, linemodel_nobroad]: - EMFit.populate_linemodel(linemodel, initial_guesses, param_bounds, log) + EMFit.populate_linemodel(linemodel_nobroad, initial_guesses, param_bounds, log) + EMFit.populate_linemodel(linemodel_broad, initial_guesses, param_bounds, log) # Initial fit - initial_linemodel_nobroad t0 = time.time() diff --git a/py/fastspecfit/fastspecfit.py b/py/fastspecfit/fastspecfit.py index 1456a69f..d690acc8 100644 --- a/py/fastspecfit/fastspecfit.py +++ b/py/fastspecfit/fastspecfit.py @@ -80,6 +80,7 @@ def parse(options=None, log=None): """ import argparse, sys + from fastspecfit.io import DEFAULT_TEMPLATEVERSION, DEFAULT_IMF parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) @@ -97,8 +98,8 @@ def parse(options=None, log=None): parser.add_argument('--constrain-age', action='store_true', help='Constrain the age of the SED.') parser.add_argument('--no-smooth-continuum', action='store_true', help='Do not fit the smooth continuum.') parser.add_argument('--percamera-models', action='store_true', help='Return the per-camera (not coadded) model spectra.') - parser.add_argument('--imf', type=str, default='chabrier', help='Initial mass function.') - parser.add_argument('--templateversion', type=str, default='1.2.0', help='Template version number.') + parser.add_argument('--imf', type=str, default=DEFAULT_IMF, help='Initial mass function.') + parser.add_argument('--templateversion', type=str, default=DEFAULT_TEMPLATEVERSION, help='Template version number.') parser.add_argument('--templates', type=str, default=None, help='Optional full path and filename to the templates.') parser.add_argument('--redrockfile-prefix', type=str, default='redrock-', help='Prefix of the input Redrock file name(s).') parser.add_argument('--specfile-prefix', type=str, default='coadd-', help='Prefix of the spectral file(s).') diff --git a/py/fastspecfit/io.py b/py/fastspecfit/io.py index 0f3454c9..4dbfb3e6 100644 --- a/py/fastspecfit/io.py +++ b/py/fastspecfit/io.py @@ -69,6 +69,8 @@ QNLINES = ['C_LYA', 'C_CIV', 'C_CIII', 'C_MgII', 'C_Hbeta', 'C_Halpha'] FLUXNORM = 1e17 # flux normalization factor for all DESI spectra [erg/s/cm2/A] +DEFAULT_TEMPLATEVERSION = '1.3.0' +DEFAULT_IMF = 'chabrier' # Taken from Redrock/0.15.4 class _ZWarningMask(object): @@ -101,8 +103,7 @@ def unpack_one_spectrum(iobj, specdata, meta, ebv, fphoto, fastphot, CTools = ContinuumTools(fphoto=fphoto, ignore_photometry=ignore_photometry) - log.info('Pre-processing object {} [targetid {} z={:.6f}].'.format( - iobj, meta[CTools.uniqueid], meta['Z'])) + log.info(f'Pre-processing object {iobj} [targetid {meta[CTools.uniqueid]} z={meta["Z"]:.6f}].') RV = 3.1 meta['EBV'] = ebv @@ -1994,14 +1995,14 @@ def select(fastfit, metadata, coadd_type, healpixels=None, tiles=None, else: return fastfit[keep], metadata[keep] -def get_templates_filename(templateversion='1.2.0', imf='chabrier'): +def get_templates_filename(templateversion=DEFAULT_TEMPLATEVERSION, imf=DEFAULT_IMF): """Get the templates filename. """ from fastspecfit.io import FTEMPLATES_DIR_NERSC templates_dir = os.path.expandvars(os.environ.get('FTEMPLATES_DIR', FTEMPLATES_DIR_NERSC)) - templates = os.path.join(templates_dir, templateversion, 'ftemplates-{}-{}.fits'.format( - imf, templateversion)) + templates = os.path.join(templates_dir, templateversion, f'ftemplates-{imf}-{templateversion}.fits') return templates + def get_qa_filename(metadata, coadd_type, outprefix=None, outdir=None, fastphot=False, log=None): """Build the QA filename. @@ -2057,8 +2058,8 @@ def _one_filename(_metadata): return pngfile -def cache_templates(templates=None, templateversion='1.2.0', imf='chabrier', - mintemplatewave=None, maxtemplatewave=40e4, vdisp_nominal=125.0, +def cache_templates(templates=None, templateversion=DEFAULT_TEMPLATEVERSION, imf=DEFAULT_IMF, + mintemplatewave=None, maxtemplatewave=40e4, vdisp_nominal=125., read_linefluxes=False, fastphot=False, log=None): """"Read the templates into a dictionary. @@ -2074,7 +2075,7 @@ def cache_templates(templates=None, templateversion='1.2.0', imf='chabrier', templates = get_templates_filename(templateversion=templateversion, imf=imf) if not os.path.isfile(templates): - errmsg = 'Templates file not found {}'.format(templates) + errmsg = f'Templates file {templates} not found.' log.critical(errmsg) raise IOError(errmsg) @@ -2083,7 +2084,7 @@ def cache_templates(templates=None, templateversion='1.2.0', imf='chabrier', templateflux = fitsio.read(templates, ext='FLUX') # [npix,nsed] templatelineflux = fitsio.read(templates, ext='LINEFLUX') # [npix,nsed] templateinfo, templatehdr = fitsio.read(templates, ext='METADATA', header=True) - + # Trim the wavelengths and select the number/ages of the templates. # https://www.sdss.org/dr14/spectro/galaxy_mpajhu if mintemplatewave is None: diff --git a/py/fastspecfit/mpi.py b/py/fastspecfit/mpi.py index 389b293f..384da21d 100644 --- a/py/fastspecfit/mpi.py +++ b/py/fastspecfit/mpi.py @@ -46,8 +46,8 @@ def get_ntargets_one(specfile, htmldir_root, outdir_root, coadd_type='healpix', def plan(comm=None, specprod=None, specprod_dir=None, coadd_type='healpix', survey=None, program=None, healpix=None, tile=None, night=None, - outdir_data='.', outdir_html='.', mp=1, merge=False, makeqa=False, - fastphot=False, overwrite=False): + outdir_data='.', mp=1, merge=False, makeqa=False, fastphot=False, + overwrite=False): import fitsio from astropy.table import Table, vstack @@ -98,7 +98,6 @@ def plan(comm=None, specprod=None, specprod_dir=None, coadd_type='healpix', outdir = os.path.join(outdir_data, specprod, subdir) htmldir = os.path.join(outdir_data, specprod, 'html', subdir) - #htmldir = os.path.join(outdir_html, specprod, subdir) def _findfiles(filedir, prefix='redrock', survey=None, program=None, healpix=None, tile=None, night=None, gzip=False): if gzip: diff --git a/py/fastspecfit/qa.py b/py/fastspecfit/qa.py index 08f08d27..3a040e37 100644 --- a/py/fastspecfit/qa.py +++ b/py/fastspecfit/qa.py @@ -1169,7 +1169,8 @@ def parse(options=None): """ import sys, argparse - + from fastspecfit.io import DEFAULT_TEMPLATEVERSION, DEFAULT_IMF + parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) parser.add_argument('--healpix', default=None, type=str, nargs='*', help="""Generate QA for all objects @@ -1201,8 +1202,8 @@ def parse(options=None): parser.add_argument('--stackfit', action='store_true', help='Generate QA for stacked spectra.') parser.add_argument('--overwrite', action='store_true', help='Overwrite existing files.') - parser.add_argument('--imf', type=str, default='chabrier', help='Initial mass function.') - parser.add_argument('--templateversion', type=str, default='1.2.0', help='Template version number.') + parser.add_argument('--imf', type=str, default=DEFAULT_IMF, help='Initial mass function.') + parser.add_argument('--templateversion', type=str, default=DEFAULT_TEMPLATEVERSION, help='Template version number.') parser.add_argument('--templates', type=str, default=None, help='Optional full path and filename to the templates.') parser.add_argument('--outprefix', default=None, type=str, help='Optional prefix for output filename.') diff --git a/py/fastspecfit/test/test_fastspecfit.py b/py/fastspecfit/test/test_fastspecfit.py index 40b21adb..af8d6a1a 100644 --- a/py/fastspecfit/test/test_fastspecfit.py +++ b/py/fastspecfit/test/test_fastspecfit.py @@ -31,8 +31,8 @@ def setUpClass(cls): cls.redrockfile = resource_filename('fastspecfit.test', 'data/redrock-4-80613-thru20210324.fits') cls.outdir = tempfile.mkdtemp() - cls.templates = os.path.join(cls.outdir, 'ftemplates-chabrier-1.2.0.fits') - cmd = 'wget -O {} https://data.desi.lbl.gov/public/external/templates/fastspecfit/1.2.0/ftemplates-chabrier-1.2.0.fits'.format(cls.templates) + cls.templates = os.path.join(cls.outdir, 'ftemplates-chabrier-1.3.0.fits') + cmd = 'wget -O {} https://data.desi.lbl.gov/public/external/templates/fastspecfit/1.3.0/ftemplates-chabrier-1.3.0.fits'.format(cls.templates) err = subprocess.call(cmd.split()) cls.cwd = os.getcwd() diff --git a/requirements.txt b/requirements.txt index 7f91679c..f035918a 100644 --- a/requirements.txt +++ b/requirements.txt @@ -7,6 +7,6 @@ seaborn matplotlib fitsio git+https://github.com/desihub/desimodel.git@0.19.0#egg=desimodel -git+https://github.com/desihub/desitarget.git@2.6.1#egg=desitarget +git+https://github.com/desihub/desitarget.git@2.7.0#egg=desitarget git+https://github.com/desihub/desispec.git@0.60.2#egg=desispec git+https://github.com/desihub/speclite.git@v0.17#egg=speclite