Skip to content

Commit

Permalink
Merge pull request #145 from desihub/vac2.1-prep
Browse files Browse the repository at this point in the history
updates and bug fixes ahead of next version of VACs
  • Loading branch information
moustakas committed Aug 7, 2023
2 parents cec13e9 + b50bd9d commit 2d2d757
Show file tree
Hide file tree
Showing 27 changed files with 1,871 additions and 2,068 deletions.
Empty file modified bin/build-fsps-templates
100644 → 100755
Empty file.
9 changes: 9 additions & 0 deletions bin/fastqa
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
#!/usr/bin/env python
"""fastspecfit QA
fastqa /global/cfs/cdirs/desi/spectro/fastspecfit/blanc/tiles/80606/deep/fastphot-0-80606-deep.fits --outdir . --ntargets 1 --firsttarget 50
"""
if __name__ == '__main__':
from fastspecfit.qa import fastqa
fastqa(comm=None)
279 changes: 3 additions & 276 deletions bin/fastspecfit-qa
Original file line number Diff line number Diff line change
@@ -1,279 +1,6 @@
#!/usr/bin/env python
"""fastspecfit QA
fastspecfit-qa --fastphotfile /global/cfs/cdirs/desi/spectro/fastspecfit/blanc/tiles/80606/deep/fastphot-0-80606-deep.fits --outdir . --ntargets 20 --firsttarget 50
"""
import pdb # for debugging
import os, sys, time
import numpy as np

from desiutil.log import get_logger
log = get_logger()

def parse(options=None):
"""Parse input arguments.
"""
import sys, argparse

parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)

parser.add_argument('--healpix', default=None, type=str, nargs='*', help="""Generate QA for all objects
with this healpixels (only defined for coadd-type 'healpix').""")
parser.add_argument('--tile', default=None, type=str, nargs='*', help='Generate QA for all objects on this tile.')
parser.add_argument('--night', default=None, type=str, nargs='*', help="""Generate QA for all objects observed on this
night (only defined for coadd-type 'pernight' and 'perexp').""")
parser.add_argument('--redux_dir', default=None, type=str, help='Top-level path to the reduced spectra.')
parser.add_argument('--redrockfiles', nargs='*', help='Optional full path to redrock file(s).')
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).')
parser.add_argument('--qnfile-prefix', type=str, default='qso_qn-', help='Prefix of the QuasarNet afterburner file(s).')
parser.add_argument('--mapdir', type=str, default=None, help='Optional directory name for the dust maps.')
parser.add_argument('--fphotodir', type=str, default=None, help='Top-level location of the source photometry.')
parser.add_argument('--fphotoinfo', type=str, default=None, help='Photometric information file.')

parser.add_argument('--emline_snrmin', type=float, default=0.0, help='Minimum emission-line S/N to be displayed.')
parser.add_argument('--nsmoothspec', type=int, default=0, help='Smoothing pixel value.')

parser.add_argument('--minspecwave', type=float, default=3500., help='Minimum spectral wavelength (Angstrom).')
parser.add_argument('--maxspecwave', type=float, default=9900., help='Maximum spectral wavelength (Angstrom).')
parser.add_argument('--minphotwave', type=float, default=0.1, help='Minimum photometric wavelength (micron).')
parser.add_argument('--maxphotwave', type=float, default=35., help='Maximum photometric wavelength (micron).')

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('--nophoto', action='store_true', help='Do not include the photometry in the model fitting.')
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.0.0', 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.')
parser.add_argument('-o', '--outdir', default='.', type=str, help='Full path to desired output directory.')

parser.add_argument('fastfitfile', nargs=1, help='Full path to fastspec or fastphot fitting results.')

if options is None:
args = parser.parse_args()
log.info(' '.join(sys.argv))
else:
args = parser.parse_args(options)
log.info('fastspecfit-qa {}'.format(' '.join(options)))

return args

def main(args=None, comm=None):
"""Main module.
"""
import fitsio
from astropy.table import Table
from fastspecfit.fastspecfit import _desiqa_one, desiqa_one
from fastspecfit.io import (DESISpectra, get_templates_filename, get_qa_filename,
read_fastspecfit, DESI_ROOT_NERSC, select)

if isinstance(args, (list, tuple, type(None))):
args = parse(args)

if args.redux_dir is None:
args.redux_dir = os.path.join(os.environ.get('DESI_ROOT', DESI_ROOT_NERSC), 'spectro', 'redux')
if not os.path.isdir(args.redux_dir):
errmsg = 'Data reduction directory {} not found.'.format(args.redux_dir)
log.critical(errmsg)
raise IOError(errmsg)

# Read the fitting results.
if not os.path.isfile(args.fastfitfile[0]):
log.warning('File {} not found.'.format(args.fastfitfile[0]))
return

fastfit, metadata, coadd_type, fastphot = read_fastspecfit(args.fastfitfile[0])

# 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))
fastfit = fastfit[keep]
metadata = metadata[keep]

fastfit, metadata = select(fastfit, metadata, coadd_type, healpixels=args.healpix,
tiles=args.tile, nights=args.night)

if coadd_type == 'custom' and args.redrockfiles is None:
errmsg = 'redrockfiles input is required if coadd_type==custom'
log.critical(errmsg)
raise IOError(errmsg)

# check for the input redshifts header card
hdr = fitsio.read_header(args.fastfitfile[0])
if 'INPUTZ' in hdr and hdr['INPUTZ']:
inputz = True
else:
inputz = False

if 'NOSCORR' in hdr and hdr['NOSCORR']:
no_smooth_continuum = True
else:
no_smooth_continuum = False

if args.outdir:
if not os.path.isdir(args.outdir):
os.makedirs(args.outdir, exist_ok=True)

pngfile = get_qa_filename(metadata, coadd_type, outprefix=args.outprefix,
outdir=args.outdir, fastphot=fastphot, log=log)

# are we overwriting?
if args.overwrite is False:
I = np.array([not os.path.isfile(_pngfile) for _pngfile in pngfile])
J = ~I
if np.sum(J) > 0:
log.info('Skipping {} existing QA files.'.format(np.sum(J)))
fastfit = fastfit[I]
metadata = metadata[I]

if len(metadata) == 0:
log.info('Done making all QA files!')
return

log.info('Building QA for {} objects.'.format(len(metadata)))

# Initialize the I/O class.

Spec = DESISpectra(redux_dir=args.redux_dir, fphotodir=args.fphotodir,
fphotoinfo=args.fphotoinfo, mapdir=args.mapdir)

templates = get_templates_filename(templateversion=args.templateversion, imf=args.imf)

def _wrap_qa(redrockfile, indx=None, stackfit=False):
if indx is None:
indx = np.arange(len(fastfit))

if stackfit:
stackids = fastfit['STACKID'][indx]
data = Spec.read_stacked(redrockfile, stackids=stackids, mp=args.mp)
args.nophoto = True # force nophoto=True

minspecwave = np.min(data[0]['coadd_wave']) - 20
maxspecwave = np.max(data[0]['coadd_wave']) + 20
else:
targetids = fastfit['TARGETID'][indx]
if inputz:
input_redshifts = fastfit['Z'][indx]
else:
input_redshifts = None
Spec.select(redrockfiles=redrockfile, targetids=targetids,
input_redshifts=input_redshifts,
redrockfile_prefix=args.redrockfile_prefix,
specfile_prefix=args.specfile_prefix,
qnfile_prefix=args.qnfile_prefix)
data = Spec.read_and_unpack(fastphot=fastphot, synthphot=True, mp=args.mp)

minspecwave = args.minspecwave
maxspecwave = args.maxspecwave

qaargs = [(data[igal], fastfit[indx[igal]], metadata[indx[igal]],
templates, coadd_type, Spec.fphoto, minspecwave, maxspecwave,
args.minphotwave, args.maxphotwave, args.emline_snrmin,
args.nsmoothspec, fastphot, args.nophoto, stackfit, inputz,
no_smooth_continuum, args.outdir, args.outprefix, log)
for igal in np.arange(len(indx))]
if args.mp > 1:
import multiprocessing
with multiprocessing.Pool(args.mp) as P:
P.map(_desiqa_one, qaargs)
else:
[desiqa_one(*_qaargs) for _qaargs in qaargs]

t0 = time.time()
if coadd_type == 'healpix':
if args.redrockfiles is not None:
for redrockfile in args.redrockfiles:
_wrap_qa(redrockfile)
else:
allspecprods = metadata['SPECPROD'].data
allsurveys = metadata['SURVEY'].data
allprograms = metadata['PROGRAM'].data
allpixels = metadata['HEALPIX'].data
for specprod in set(allspecprods):
for survey in set(allsurveys):
for program in set(allprograms):
for pixel in set(allpixels):
indx = np.where((specprod == allspecprods) * (survey == allsurveys) *
(program == allprograms) * (pixel == allpixels))[0]
if len(indx) == 0:
#log.warning('No object found with specprod={}, survey={}, program={}, and healpixel={}!'.format(
# specprod, survey, program, pixel))
continue
redrockfile = os.path.join(args.redux_dir, specprod, 'healpix', str(survey), str(program), str(pixel // 100),
str(pixel), 'redrock-{}-{}-{}.fits'.format(survey, program, pixel))
_wrap_qa(redrockfile, indx)
elif coadd_type == 'custom':
for redrockfile in args.redrockfiles:
_wrap_qa(redrockfile)
elif coadd_type == 'stacked':
for redrockfile in args.redrockfiles:
_wrap_qa(redrockfile, stackfit=True)
else:
allspecprods = metadata['SPECPROD'].data
alltiles = metadata['TILEID'].astype(str).data
allnights = metadata['NIGHT'].astype(str).data
allpetals = metadata['FIBER'].data // 500
if coadd_type == 'cumulative':
for specprod in set(allspecprods):
for tile in set(alltiles):
for petal in set(allpetals):
indx = np.where((specprod == allspecprods) * (tile == alltiles) * (petal == allpetals))[0]
if len(indx) == 0:
#log.warning('No object found with tileid={} and petal={}!'.format(
# tile, petal))
continue
redrockfile = os.path.join(args.redux_dir, specprod, 'tiles', 'cumulative', str(tile), allnights[indx[0]],
'redrock-{}-{}-thru{}.fits'.format(petal, tile, allnights[indx[0]]))
_wrap_qa(redrockfile, indx)
elif coadd_type == 'pernight':
for specprod in set(allspecprods):
for night in set(allnights):
for tile in set(alltiles):
for petal in set(allpetals):
indx = np.where((specprod == allspecprods) * (night == allnights) *
(tile == alltiles) * (petal == allpetals))[0]
if len(indx) == 0:
continue
redrockfile = os.path.join(args.redux_dir, specprod, 'tiles', 'pernight', str(tile), str(night),
'redrock-{}-{}-{}.fits'.format(petal, tile, night))
_wrap_qa(redrockfile, indx)
elif coadd_type == 'perexp':
allexpids = metadata['EXPID'].data
for specprod in set(allspecprods):
for night in set(allnights):
for expid in set(allexpids):
for tile in set(alltiles):
for petal in set(allpetals):
indx = np.where((specprod == allspecprods) * (night == allnights) *
(expid == allexpids) * (tile == alltiles) *
(petal == allpetals))[0]
if len(indx) == 0:
continue
redrockfile = os.path.join(args.redux_dir, specprod, 'tiles', 'perexp', str(tile), '{:08d}'.format(expid),
'redrock-{}-{}-exp{:08d}.fits'.format(petal, tile, expid))
_wrap_qa(redrockfile, indx)

log.info('QA for everything took: {:.2f} sec'.format(time.time()-t0))

if __name__ == '__main__':
main()

from desiutil.log import get_logger
log = get_logger()
log.warning('fastspecfit-qa is deprecated; please use fastqa')
2 changes: 1 addition & 1 deletion bin/fastspecfit-setup.sh
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ elif [[ $1 == "env" ]]; then

export DESI_ROOT='/global/cfs/cdirs/desi'
export DUST_DIR='/global/cfs/cdirs/cosmo/data/dust/v0_1'
export DR9_DIR='/global/cfs/cdirs/desi/external/legacysurvey/dr9'
export FPHOTO_DIR='/global/cfs/cdirs/desi/external/legacysurvey/dr9'
export FTEMPLATES_DIR='/global/cfs/cdirs/desi/external/templates/fastspecfit'

export OMP_NUM_THREADS=1
Expand Down
21 changes: 10 additions & 11 deletions bin/get-cutouts
Original file line number Diff line number Diff line change
Expand Up @@ -135,10 +135,9 @@ def plan(comm=None, specprod=None, coadd_type='healpix',
else:
rank, size = comm.rank, comm.size

desi_root = '/global/cfs/cdirs/desi'
#desi_root = os.environ.get('DESI_ROOT', DESI_ROOT_NERSC)
# look for data in the standard location
desi_root = os.environ.get('DESI_ROOT', '/dvs_ro/cfs/cdirs/desi')

# look for data in the standard location
if coadd_type == 'healpix':
subdir = 'healpix'
if healpix is not None:
Expand Down Expand Up @@ -352,7 +351,7 @@ def main():
parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser.add_argument('--coadd-type', type=str, default='healpix', choices=['healpix', 'cumulative', 'pernight', 'perexp'],
help='Specify which type of spectra/zbest files to process.')
parser.add_argument('--specprod', type=str, default='fuji', #choices=['fuji', 'guadalupe', 'iron'],
parser.add_argument('--specprod', type=str, default='iron', #choices=['fuji', 'guadalupe', 'iron'],
help='Spectroscopic production to process.')

parser.add_argument('--healpix', type=str, default=None, help='Comma-separated list of healpixels to process.')
Expand All @@ -367,13 +366,14 @@ def main():
parser.add_argument('--nompi', action='store_true', help='Do not use MPI parallelism.')
parser.add_argument('--dry-run', action='store_true', help='Generate but do not run commands.')

parser.add_argument('--outdir-html', default='/global/cfs/cdirs/desi/users/ioannis/fastspecfit',
type=str, help='Base output HTML directory.')
parser.add_argument('--outdir-data', default='/global/cfs/cdirs/desi/spectro/fastspecfit',
type=str, help='Base output data directory.')
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.')

args = parser.parse_args()

outdir_data = os.path.expandvars(args.outdir_data)
outdir_html = os.path.expandvars(args.outdir_html)

if args.nompi:
comm = None
else:
Expand All @@ -397,11 +397,10 @@ def main():
plan(comm=comm, specprod=args.specprod, coadd_type=args.coadd_type,
survey=args.survey, program=args.program,
healpix=args.healpix, tile=args.tile, night=args.night,
outdir_data=args.outdir_data, outdir_html=args.outdir_html,
outdir_data=outdir_data, outdir_html=outdir_html,
overwrite=args.overwrite, mp=args.mp)
else:
do_cutouts(args, comm=comm, outdir_data=args.outdir_data,
outdir_html=args.outdir_html)
do_cutouts(args, comm=comm, outdir_data=outdir_data, outdir_html=outdir_html)

if __name__ == '__main__':
main()

0 comments on commit 2d2d757

Please sign in to comment.