diff --git a/bin/wrap-fastframe b/bin/wrap-fastframe index 261f1e665..dbe9bc617 100755 --- a/bin/wrap-fastframe +++ b/bin/wrap-fastframe @@ -28,7 +28,7 @@ else: rank = 0 #- The proceed with other imports -import sys, os, glob, time +import sys, os, glob, time, subprocess tstart = time.time() import numpy as np @@ -89,8 +89,24 @@ if rank < len(simspecfiles): try: t0 = time.time() - desisim.scripts.fastframe.main(cmd.split()[1:]) + logfile = filename.replace('simspec', 'fastframe').replace('.fits', '.log') + assert logfile != filename + + #- Use subprocess.call instead of fastframe.main() to avoid + #- potential memory leaks and separate logging; this does incur + #- extra python interpreter startup time. + ### desisim.scripts.fastframe.main(cmd.split()[1:]) + + print('logging to {}'.format(logfile)) + with open(logfile, 'w') as logx: + err = subprocess.call(cmd.split(), stdout=logx, stderr=logx) + runtime = time.time() - t0 + if err != 0: + print("ERROR: rank {} simspec {} error code {} after {:.1f} sec".format( + rank, os.path.basename(filename), err, runtime)) + raise RuntimeError + print("rank {} took {:.1f} seconds for {} frame".format(rank, runtime, flavor)) sys.stdout.flush() except: diff --git a/bin/wrap-newexp b/bin/wrap-newexp index 191b560c2..d23ae98b9 100755 --- a/bin/wrap-newexp +++ b/bin/wrap-newexp @@ -8,6 +8,10 @@ MPI wrapper for newexp #- from login node where MPI doesn't work from __future__ import absolute_import, division, print_function import argparse +import subprocess +import sys + +import desisim.io parser=argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) @@ -49,7 +53,7 @@ from astropy.table import Table from desisim.util import dateobs2night import desisim.io -import desisim.scripts.newexp +import desisim.scripts.newexp_mock import desisim.scripts.newflat import desisim.scripts.newarc @@ -166,7 +170,7 @@ if rank < len(todo): (flavor, night, expid, obsnum, obscond) = thisobs assert flavor in ('science', 'arc', 'flat') if flavor == 'science': - cmd = "newexp --obslist {} --fiberassign {} --mockdir {}".format( + cmd = "newexp-mock --obslist {} --fiberassign {} --mockdir {}".format( args.obslist, args.fiberassign, args.mockdir) if args.outdir is not None: cmd = cmd + " --outdir {}".format(args.outdir) @@ -183,19 +187,33 @@ if rank < len(todo): #- TODO: what about logging? Use desisim.pipe log redirection logic? print('RUNNING: {}'.format(cmd)) + sys.stdout.flush() if args.dryrun: continue try: t0 = time.time() - if cmd.startswith('newexp'): - desisim.scripts.newexp.main(cmd.split()[1:]) - elif cmd.startswith('newflat'): - desisim.scripts.newflat.main(cmd.split()[1:]) - elif cmd.startswith('newarc'): - desisim.scripts.newarc.main(cmd.split()[1:]) - else: - print('ERROR: unknown command {}'.format(cmd)) + logfile = desisim.io.findfile( + 'simspec', night, expid).replace('.fits', '.log') + print('logging {}-{} {} to {}'.format(night, expid, flavor, logfile)) + + #- Spawn call to avoid memory leak problems from repeated + #- desisim.scripts.newexp_mock.main() calls within same interpreter + with open(logfile, 'w') as logx: + err = subprocess.call(cmd.split(), stdout=logx, stderr=logx) + # if cmd.startswith('newexp'): + # desisim.scripts.newexp_mock.main(cmd.split()[1:]) + # elif cmd.startswith('newflat'): + # desisim.scripts.newflat.main(cmd.split()[1:]) + # elif cmd.startswith('newarc'): + # desisim.scripts.newarc.main(cmd.split()[1:]) + # else: + # print('ERROR: unknown command {}'.format(cmd)) + if err != 0: + print('ERROR: night {} expid {} failed with error {}'.format( + night, expid, err)) + raise RuntimeError + runtime = time.time() - t0 print("{} took {:.1f} seconds".format(cmd.split()[0], runtime)) sys.stdout.flush() diff --git a/doc/changes.rst b/doc/changes.rst index 7bc532a6c..f2e174f50 100644 --- a/doc/changes.rst +++ b/doc/changes.rst @@ -18,11 +18,13 @@ desisim change log * Miscellaneous polishing in QA (velocity, clip before RMS, extend [OII] flux, S/N per Ang) * Bug fix: correctly select both "bright" and "faint" BGS templates by default (`PR #257`_). +* Updates for newexp/fastframe wrappers for end-to-end sims (`PR #258`_). .. _`PR #250`: https://github.com/desihub/desisim/pull/250 .. _`PR #252`: https://github.com/desihub/desisim/pull/252 .. _`PR #254`: https://github.com/desihub/desisim/pull/254 .. _`PR #257`: https://github.com/desihub/desisim/pull/257 +.. _`PR #258`: https://github.com/desihub/desisim/pull/258 0.20.0 (2017-07-12) ------------------- diff --git a/py/desisim/io.py b/py/desisim/io.py index fa27e5c3f..38215adfd 100644 --- a/py/desisim/io.py +++ b/py/desisim/io.py @@ -60,6 +60,8 @@ def findfile(filetype, night, expid, camera=None, outdir=None, mkdir=True): simpix = '{outdir:s}/simpix-{expid:08d}.fits', simfibermap = '{outdir:s}/fibermap-{expid:08d}.fits', pix = '{outdir:s}/pix-{camera:s}-{expid:08d}.fits', + fastframelog = '{outdir:s}/fastframe-{expid:08d}.log', + newexplog = '{outdir:s}/newexp-{expid:08d}.log', ) #- Do we know about this kind of file? @@ -129,9 +131,15 @@ def write_simspec(sim, truth, fibermap, obs, expid, night, outdir=None, filename header['NIGHT'] = night header['EXPTIME'] = sim.observation.exposure_time.to('s').value if obs is not None: - for key in obs.keys(): - if key not in header: - header[key] = obs[key] + try: + keys = obs.keys() + except AttributeError: + keys = obs.dtype.names + + for key in keys: + shortkey = key[0:8] #- FITS keywords can only be 8 char + if shortkey not in header: + header[shortkey] = obs[key] if 'DOSVER' not in header: header['DOSVER'] = 'SIM' if 'FEEVER' not in header: @@ -880,7 +888,7 @@ def simdir(night='', mkdir=False): Return $DESI_SPECTRO_SIM/$PIXPROD/{night} If mkdir is True, create directory if needed """ - dirname = os.path.join(os.getenv('DESI_SPECTRO_SIM'), os.getenv('PIXPROD'), night) + dirname = os.path.join(os.getenv('DESI_SPECTRO_SIM'), os.getenv('PIXPROD'), str(night)) if mkdir and not os.path.exists(dirname): os.makedirs(dirname) diff --git a/py/desisim/obs.py b/py/desisim/obs.py index 3b14b5e0c..c5fa56e1e 100644 --- a/py/desisim/obs.py +++ b/py/desisim/obs.py @@ -175,7 +175,7 @@ def new_exposure(program, nspec=5000, night=None, expid=None, tileid=None, if exptime is not None: obsconditions['EXPTIME'] = exptime - sim = simulate_spectra(wave, 1e-17*flux, fibermap=fibermap, obsconditions=obsconditions) + sim = simulate_spectra(wave, flux, fibermap=fibermap, obsconditions=obsconditions) #- Override $DESI_SPECTRO_DATA in order to write to simulation area datadir_orig = os.getenv('DESI_SPECTRO_DATA') @@ -444,7 +444,9 @@ def update_obslog(obstype='science', program='DARK', expid=None, dateobs=None, INSERT OR REPLACE INTO obslog(expid,dateobs,night,obstype,program,tileid,ra,dec) VALUES (?,?,?,?,?,?,?,?) """ - db.execute(insert, (expid, time.mktime(dateobs), night, obstype.upper(), program.upper(), tileid, ra, dec)) + db.execute(insert, (int(expid), time.mktime(dateobs), str(night), + str(obstype.upper()), str(program.upper()), int(tileid), + float(ra), float(dec))) db.commit() return expid, dateobs diff --git a/py/desisim/scripts/fastframe.py b/py/desisim/scripts/fastframe.py index c57a04b5b..1dc4b9283 100644 --- a/py/desisim/scripts/fastframe.py +++ b/py/desisim/scripts/fastframe.py @@ -22,8 +22,10 @@ def parse(options=None): parser = argparse.ArgumentParser(usage = "{prog} [options]") parser.add_argument("--simspec", type=str, help="input simspec file") parser.add_argument("--outdir", type=str, help="output directory") - parser.add_argument("--firstspec", type=int, default=0, help="first spectrum to simulate") - parser.add_argument("--nspec", type=int, default=5000, help="number of spectra to simulate") + parser.add_argument("--firstspec", type=int, default=0, + help="first spectrum to simulate") + parser.add_argument("--nspec", type=int, default=5000, + help="number of spectra to simulate") if options is None: args = parser.parse_args() @@ -61,12 +63,13 @@ def main(args=None): flux = simspec.flux ii = slice(firstspec, firstspec+nspec) if simspec.flavor == 'science': - sim = desisim.simexp.simulate_spectra(wave, 1e-17*flux[ii], fibermap=fibermap[ii], - obsconditions=obs, dwave_out=1.0) + sim = desisim.simexp.simulate_spectra(wave, flux[ii], + fibermap=fibermap[ii], obsconditions=obs, dwave_out=1.0) elif simspec.flavor in ['arc', 'flat', 'calib']: x = fibermap['X_TARGET'] y = fibermap['Y_TARGET'] - fiber_area = desisim.simexp.fiber_area_arcsec2(fibermap['X_TARGET'], fibermap['Y_TARGET']) + fiber_area = desisim.simexp.fiber_area_arcsec2( + fibermap['X_TARGET'], fibermap['Y_TARGET']) surface_brightness = (flux.T / fiber_area).T config = desisim.simexp._specsim_config_for_wave(wave, dwave_out=1.0) # sim = specsim.simulator.Simulator(config, num_fibers=nspec) @@ -74,7 +77,8 @@ def main(args=None): sim.observation.exposure_time = simspec.header['EXPTIME'] * u.s sbunit = 1e-17 * u.erg / (u.Angstrom * u.s * u.cm ** 2 * u.arcsec ** 2) xy = np.vstack([x, y]).T * u.mm - sim.simulate(calibration_surface_brightness=surface_brightness[ii]*sbunit, focal_positions=xy[ii]) + sim.simulate(calibration_surface_brightness=surface_brightness[ii]*sbunit, + focal_positions=xy[ii]) else: raise ValueError('Unknown simspec flavor {}'.format(simspec.flavor)) @@ -89,7 +93,8 @@ def main(args=None): results['random_noise_electrons']).T ivar = 1.0 / results['variance_electrons'].T R = Resolution(sim.instrument.cameras[i].get_output_resolution_matrix()) - Rdata = np.tile(R.data.T, nspec).T.reshape(nspec, R.data.shape[0], R.data.shape[1]) + Rdata = np.tile(R.data.T, nspec).T.reshape( + nspec, R.data.shape[0], R.data.shape[1]) assert np.all(Rdata[0] == R.data) assert phot.shape == (nspec, len(wave)) for spectro in range(10): @@ -105,6 +110,7 @@ def main(args=None): meta['CAMERA'] = camera frame = Frame(wave, xphot, xivar, resolution_data=Rdata[0:imax-imin], spectrograph=spectro, fibermap=xfibermap, meta=meta) - outfile = desispec.io.findfile('frame', night, expid, camera, outdir=args.outdir) + outfile = desispec.io.findfile('frame', night, expid, camera, + outdir=args.outdir) print('writing {}'.format(outfile)) desispec.io.write_frame(outfile, frame) diff --git a/py/desisim/scripts/newexp_mock.py b/py/desisim/scripts/newexp_mock.py index e1cbab878..fb5ba117e 100644 --- a/py/desisim/scripts/newexp_mock.py +++ b/py/desisim/scripts/newexp_mock.py @@ -101,7 +101,12 @@ def main(args=None): args.nspec = len(fiberassign) log.info('Simulating night {} expid {} tile {}'.format(night, args.expid, tileid)) - flux, wave, meta = get_mock_spectra(fiberassign, mockdir=args.mockdir) + try: + flux, wave, meta = get_mock_spectra(fiberassign, mockdir=args.mockdir) + except Exception as err: + log.fatal('Failed obsnum {} fiberassign {} tile {}'.format(args.obsnum, args.fiberassign, tileid)) + raise err + sim, fibermap = simscience((flux, wave, meta), fiberassign, obsconditions=obs, nspec=args.nspec) #- TODO: header keyword code is replicated from obs.new_exposure() @@ -114,17 +119,17 @@ def main(args=None): FLAVOR = ('science', 'Flavor [arc, flat, science, zero, ...]'), TELRA = (telera, 'Telescope pointing RA [degrees]'), TELDEC = (teledec, 'Telescope pointing dec [degrees]'), - AIRMASS = (obsconditions['AIRMASS'], 'Airmass at middle of exposure'), - EXPTIME = (obsconditions['EXPTIME'], 'Exposure time [sec]'), - SEEING = (obsconditions['SEEING'], 'Seeing FWHM [arcsec]'), - MOONFRAC = (obsconditions['MOONFRAC'], 'Moon illumination fraction 0-1; 1=full'), - MOONALT = (obsconditions['MOONALT'], 'Moon altitude [degrees]'), - MOONSEP = (obsconditions['MOONSEP'], 'Moon:tile separation angle [degrees]'), + AIRMASS = (obs['AIRMASS'], 'Airmass at middle of exposure'), + EXPTIME = (obs['EXPTIME'], 'Exposure time [sec]'), + SEEING = (obs['SEEING'], 'Seeing FWHM [arcsec]'), + MOONFRAC = (obs['MOONFRAC'], 'Moon illumination fraction 0-1; 1=full'), + MOONALT = (obs['MOONALT'], 'Moon altitude [degrees]'), + MOONSEP = (obs['MOONSEP'], 'Moon:tile separation angle [degrees]'), ) header['DATE-OBS'] = (sim.observation.exposure_start.isot, 'Start of exposure') #- Write fibermap to $DESI_SPECTRO_SIM/$PIXPROD not $DESI_SPECTRO_DATA - fibermap.meta.extend(header) + fibermap.meta.update(header) fibermap.write(desisim.io.findfile('simfibermap', night, args.expid, outdir=args.outdir), overwrite=args.clobber) diff --git a/py/desisim/simexp.py b/py/desisim/simexp.py index f72f509ba..4ea4a2f9d 100644 --- a/py/desisim/simexp.py +++ b/py/desisim/simexp.py @@ -10,6 +10,7 @@ import astropy.table import astropy.time from astropy.io import fits +import fitsio import desitarget import desispec.io @@ -290,18 +291,18 @@ def fibermeta2fibermap(fiberassign, meta): fibermap[c] = fiberassign[c] #- MAG and FILTER; ignore warnings from negative flux - #- these are deprecated anyway and will be replaced with FLUX_G, FLUX_R, etc. + #- these are deprecated anyway and will be replaced with FLUX_G, FLUX_R, + #- etc. in the fibermap as well with warnings.catch_warnings(): warnings.simplefilter('ignore') - ugrizy = 22.5 - 2.5 * np.log10(meta['DECAM_FLUX'].data) - wise = 22.5 - 2.5 * np.log10(meta['WISE_FLUX'].data) - fibermap['FILTER'][:, :5] = ['DECAM_G', 'DECAM_R', 'DECAM_Z', 'WISE_W1', 'WISE_W2'] - fibermap['MAG'][:, 0] = ugrizy[:, 1] - fibermap['MAG'][:, 1] = ugrizy[:, 2] - fibermap['MAG'][:, 2] = ugrizy[:, 4] - fibermap['MAG'][:, 3] = wise[:, 0] - fibermap['MAG'][:, 4] = wise[:, 1] + fibermap['FILTER'][:, :5] = \ + ['DECAM_G', 'DECAM_R', 'DECAM_Z', 'WISE_W1', 'WISE_W2'] + fibermap['MAG'][:, 0] = 22.5 - 2.5 * np.log10(meta['FLUX_G'].data) + fibermap['MAG'][:, 1] = 22.5 - 2.5 * np.log10(meta['FLUX_R'].data) + fibermap['MAG'][:, 2] = 22.5 - 2.5 * np.log10(meta['FLUX_Z'].data) + fibermap['MAG'][:, 3] = 22.5 - 2.5 * np.log10(meta['FLUX_W1'].data) + fibermap['MAG'][:, 4] = 22.5 - 2.5 * np.log10(meta['FLUX_W2'].data) #- set OBJTYPE #- TODO: what about MWS science targets that are also standard stars? @@ -338,7 +339,8 @@ def simulate_spectra(wave, flux, fibermap=None, obsconditions=None, dwave_out=No Args: wave (array): 1D wavelengths in Angstroms - flux (array): 2D[nspec,nwave] flux in erg/s/cm2/Angstrom + flux (array): 2D[nspec,nwave] flux in 1e-17 erg/s/cm2/Angstrom + or astropy Quantity with flux units Optional: fibermap: table from fiberassign or fibermap; uses X/YFOCAL_DESIGN, TARGETID, DESI_TARGET @@ -362,7 +364,7 @@ def simulate_spectra(wave, flux, fibermap=None, obsconditions=None, dwave_out=No #- Convert to unit-ful quantities for specsim if not isinstance(flux, u.Quantity): - fluxunits = u.erg / (u.Angstrom * u.s * u.cm**2) + fluxunits = 1e-17 * u.erg / (u.Angstrom * u.s * u.cm**2) flux = flux * fluxunits if not isinstance(wave, u.Quantity): @@ -635,37 +637,59 @@ def read_mock_spectra(truthfile, targetids, mockdir=None): wave[nwave]: wavelengths in Angstroms truth[nspec]: metadata truth table ''' - with fits.open(truthfile, memmap=False) as fx: - truth = fx['TRUTH'].data - wave = fx['WAVE'].data - flux = fx['FLUX'].data + if len(targetids) != len(np.unique(targetids)): + from desiutil.log import get_logger + log = get_logger() + log.error("Requested TARGETIDs for {} are not unique".format( + os.path.basename(truthfile))) + + #- astropy.io.fits doesn't return a real ndarray, causing problems + #- with the reordering downstream so use fitsio instead + # with fits.open(truthfile, memmap=False) as fx: + # truth = fx['TRUTH'].data + # wave = fx['WAVE'].data + # flux = fx['FLUX'].data + with fitsio.FITS(truthfile) as fx: + truth = fx['TRUTH'].read() + wave = fx['WAVE'].read() + flux = fx['FLUX'].read() + + missing = np.in1d(targetids, truth['TARGETID'], invert=True) + if np.any(missing): + missingids = targetids[missing] + raise ValueError('Targets missing from {}: {}'.format(truthfile, missingids)) #- Trim to just the spectra for these targetids ii = np.in1d(truth['TARGETID'], targetids) - if np.count_nonzero(ii) != len(targetids): - jj = ~np.in1d(targetids, truth['TARGETID']) - missingids = targetids[jj] - raise ValueError('Targets missing from {}: {}'.format(truthfile, missingids)) - flux = flux[ii] truth = truth[ii] assert set(targetids) == set(truth['TARGETID']) #- sort truth to match order of input targetids - #- Sort tiles to match order of tileids - i = np.argsort(targetids) - j = np.argsort(truth['TARGETID']) - k = np.argsort(i) - - truth = truth[j[k]] - flux = flux[j[k]] - assert np.all(truth['TARGETID'] == targetids) + if len(targetids) == len(truth['TARGETID']): + i = np.argsort(targetids) + j = np.argsort(truth['TARGETID']) + k = np.argsort(i) + reordered_truth = truth[j[k]] + reordered_flux = flux[j[k]] + else: + #- Slower, but works even with repeated TARGETIDs + ii = np.argsort(truth['TARGETID']) + sorted_truthids = truth['TARGETID'][ii] + reordered_flux = np.empty(shape=(len(targetids), flux.shape[1]), dtype=flux.dtype) + reordered_truth = np.empty(shape=(len(targetids),), dtype=truth.dtype) + for j, tx in enumerate(targetids): + k = np.searchsorted(sorted_truthids, tx) + reordered_flux[j] = flux[ii[k]] + reordered_truth[j] = truth[ii[k]] + + assert np.all(reordered_truth['TARGETID'] == targetids) wave = desispec.io.util.native_endian(wave).astype(np.float64) - flux = desispec.io.util.native_endian(flux).astype(np.float64) + reordered_flux = desispec.io.util.native_endian(reordered_flux).astype(np.float64) - return flux, wave, truth + return reordered_flux, wave, reordered_truth def targets2truthfiles(targets, basedir, nside=64): ''' @@ -699,7 +723,7 @@ def targets2truthfiles(targets, basedir, nside=64): filename = mockio.findfile('truth', nside, ipix, basedir=basedir) truthfiles.append(filename) ii = (pixels == ipix) - targetids.append(targets['TARGETID'][ii]) + targetids.append(np.asarray(targets['TARGETID'][ii])) return truthfiles, targetids