changes from master

.travis.yml

@@ -45,16 +45,15 @@ env:
         # to repeat them for all configurations.
         # - NUMPY_VERSION=1.10
         # - SCIPY_VERSION=0.16
-        - ASTROPY_VERSION=2.0.4
+        - ASTROPY_VERSION=2.0.7
         # - SPHINX_VERSION=1.6.6
-        - DESIUTIL_VERSION=1.9.9
+        - DESIUTIL_VERSION=1.9.11
         - SPECLITE_VERSION=0.7
         - SPECTER_VERSION=0.8.3
-        - DESIMODEL_VERSION=0.9.2
-        # - DESIMODEL_DATA=branches/test-0.9
-        - DESIMODEL_DATA=tags/0.9.2
+        - DESIMODEL_VERSION=0.9.6
+        - DESIMODEL_DATA=branches/test-0.9.6
         - DESISIM_TESTDATA_VERSION=master
-        - SPECSIM_VERSION=v0.11
+        - SPECSIM_VERSION=v0.12
         # - DESISPEC_VERSION=0.18.0
         - DESISPEC_VERSION=master
         - DESITARGET_VERSION=0.19.0

bin/wrap-fastframe

@@ -54,7 +54,7 @@ else:
 
 #- Assemble full list of simspec files; group by flavor to help MPI balancing
 if rank == 0:
-    allfiles = sorted(glob.glob('{}/*/simspec*.fits'.format(desisim.io.simdir())))
+    allfiles = sorted(glob.glob('{}/*/*/simspec*.fits'.format(desisim.io.simdir())))
 
     #- Generate list of cameras
     cameras = list()
@@ -64,7 +64,7 @@ if rank == 0:
 
     simspecfiles = list()
     for filename in allfiles:
-        night = os.path.basename(os.path.split(filename)[0])
+        night = os.path.basename(os.path.dirname(os.path.dirname(filename)))
         if (args.start <= night) & (night < args.stop):
             flavor = fits.getval(filename, 'FLAVOR')
             expid = fits.getval(filename, 'EXPID')

doc/changes.rst

@@ -2,7 +2,24 @@
 desisim change log
 ==================
 
-0.28.0 (unreleased)
+0.29.1 (unreleased)
+-------------------
+
+* No changes yet.
+
+0.29.0 (2018-07-26)
+-------------------
+
+* Option in quickspectra to write the full sim table (`PR #392`_).
+* Option to use Gaussian instead of Poisson for QSO DLA.
+  Requires specsim >= v0.12 (`PR #393`_).
+* Use `overwrite` instead of `clobber` for `astropy.io.fits` (`PR #395`_).
+
+.. _`PR #392`: https://github.com/desihub/desisim/pull/392
+.. _`PR #393`: https://github.com/desihub/desisim/pull/393
+.. _`PR #395`: https://github.com/desihub/desisim/pull/395
+
+0.28.0 (2018-07-18)
 -------------------
 
 * Add BALs to templates.QSO class (`PR #321`_).
@@ -20,6 +37,9 @@ desisim change log
 * Add rest-frame option to templates.SIMQSO (`PR #377`_).
 * Optionally change output wave vector in templates.SIMQSO when noresample=True
   or restframe=True (`PR #383`_).
+* Fix ``newexp-mock`` and ``wrap-fastframe`` file parsing for ``NIGHT/EXPID/*.*``
+  vs. ``NIGHT/*.*``.
+* Speed up emission line simulation when using ``MKL >= 2018.0.2`` (`PR #390`_).
 
 .. _`PR #321`: https://github.com/desihub/desisim/pull/321
 .. _`PR #349`: https://github.com/desihub/desisim/pull/349
@@ -35,6 +55,7 @@ desisim change log
 .. _`PR #370`: https://github.com/desihub/desisim/pull/370
 .. _`PR #377`: https://github.com/desihub/desisim/pull/377
 .. _`PR #383`: https://github.com/desihub/desisim/pull/383
+.. _`PR #390`: https://github.com/desihub/desisim/pull/390
 
 0.27.0 (2018-03-29)
 -------------------
@@ -173,7 +194,7 @@ Requires desitarget < 0.19.0
 * Add BGS, MWS to z_find QA
 * 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`_).  
+  (`PR #257`_).
 * Updates for newexp/fastframe wrappers for end-to-end sims (`PR #258`_).
 
 .. _`PR #250`: https://github.com/desihub/desisim/pull/250

py/desisim/_version.py

@@ -1 +1 @@
-__version__ = '0.27.0.dev1392'
+__version__ = '0.29.0.dev1617'

py/desisim/dla.py

@@ -37,7 +37,7 @@ def insert_dlas(wave, zem, rstate=None, seed=None, fNHI=None, debug=False, **kwa
     log = get_logger()
     # Init
     if rstate is None:
-        rstate = np.random.RandomState(seed)
+        rstate = np.random.RandomState(seed) # this is breaking the chain of randoms if seed is None
     if fNHI is None:
         fNHI = init_fNHI(**kwargs)
 

py/desisim/io.py

@@ -226,7 +226,7 @@ def write_simspec(sim, truth, fibermap, obs, expid, night, outdir=None, filename
         hx.append(obs_hdu)
 
     log.info('Writing {}'.format(filename))
-    hx.writeto(filename, clobber=overwrite)
+    hx.writeto(filename, overwrite=overwrite)
 
 def write_simspec_arc(filename, wave, phot, header, fibermap, overwrite=False):
     '''
@@ -263,7 +263,7 @@ def write_simspec_arc(filename, wave, phot, header, fibermap, overwrite=False):
     hx.append(fibermap_hdu)
 
     log.info('Writing {}'.format(filename))
-    hx.writeto(filename, clobber=overwrite)
+    hx.writeto(filename, overwrite=overwrite)
     return filename
 
 class SimSpecCamera(object):
@@ -992,21 +992,29 @@ def write_templates(outfile, flux, wave, meta):
     hx.append(hdu_meta)
 
     log.info('Writing {}'.format(outfile))
-    try:
-        hx.writeto(outfile, overwrite=True)
-    except:
-        hx.writeto(outfile, clobber=True)
+    hx.writeto(outfile, overwrite=True)
 
 #-------------------------------------------------------------------------
 #- Utility functions
 
-def simdir(night='', expid=0, mkdir=False):
+def simdir(night=None, expid=None, 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'),
-         str(night), '{:08d}'.format(expid))
+
+    if (night is None) and (expid is None):
+        dirname = os.path.join(
+            os.getenv('DESI_SPECTRO_SIM'), os.getenv('PIXPROD')
+            )
+    #- must provide night+expid, and catch old usage where mkdir was 2nd arg
+    elif (expid is None) or isinstance(expid, bool):
+        raise ValueError("Must provide int expid, not {}".format(expid))
+    else:
+        dirname = os.path.join(
+            os.getenv('DESI_SPECTRO_SIM'), os.getenv('PIXPROD'),
+            str(night), '{:08d}'.format(expid)
+            )
     if mkdir and not os.path.exists(dirname):
         os.makedirs(dirname)
 

py/desisim/qso_template/desi_qso_templ.py

@@ -419,7 +419,7 @@ def desi_qso_templates(z_wind=0.2, zmnx=(0.4,4.), outfil=None, N_perz=500,
     tbhdu.header.set('EXTNAME','METADATA')
 
     hdulist = fits.HDUList([hdu, tbhdu])
-    hdulist.writeto(outfil, clobber=True)
+    hdulist.writeto(outfil, overwrite=True)
 
     return final_wave, final_spec, final_z
 
@@ -546,7 +546,7 @@ def repackage_coeff(boss_pca_fil=None, sdss_pca_fil=None,
 
     hdulist = fits.HDUList([phdu, bp_hdu, bz_hdu, sp_hdu, sz_hdu,
                             e_hdu, ew_hdu])
-    hdulist.writeto(outfil, clobber=True)
+    hdulist.writeto(outfil, overwrite=True)
     print('Wrote {:s}'.format(outfil))
 
 

py/desisim/qso_template/fit_boss_qsos.py

@@ -322,7 +322,7 @@ def failed_parallel():
             outfil = 'BOSS_DR10Lya_PCA_values_nocut.fits'
         else:
             outfil = 'BOSS_DR10Lya_PCA_values.fits'
-    thdulist.writeto(outfil, clobber=True)
+    thdulist.writeto(outfil, overwrite=True)
 
     # Done
     #xdb.set_trace()
@@ -372,7 +372,7 @@ def splice_fits(flg=0):
     table_hdu = fits.BinTableHDU.from_columns(np.array(full_tab.filled()))
     thdulist = fits.HDUList([prihdu, table_hdu])
     print('Writing {:s} table, with {:d} rows'.format(outfil,len(full_tab)))
-    thdulist.writeto(outfil, clobber=True)
+    thdulist.writeto(outfil, overwrite=True)
 
 
 ## ################
@@ -415,7 +415,7 @@ def splice_fits(flg=0):
     prihdu = fits.PrimaryHDU(header=prihdr)
 
     thdulist = fits.HDUList([prihdu, tbhdu])
-    thdulist.writeto(outfil, clobber=True)
+    thdulist.writeto(outfil, overwrite=True)
 
     # Done
     #xdb.set_trace()

py/desisim/scripts/newexp_mock.py

@@ -59,7 +59,7 @@ def main(args=None):
     fiberassign = astropy.table.Table.read(args.fiberassign, 'FIBERASSIGN')
 
     if args.outdir is None:
-        args.outdir = desisim.io.simdir(night=night, mkdir=True)
+        args.outdir = desisim.io.simdir(night=night, expid=args.expid, mkdir=True)
 
     if args.nspec is None:
         args.nspec = len(fiberassign)

py/desisim/scripts/quickquasars.py

@@ -208,9 +208,9 @@ def simulate_one_healpix(ifilename,args,model,obsconditions,decam_and_wise_filte
 ##ALMA:added to set transmission to 1 for z>zqso, this can be removed when transmission is corrected.        
     for ii in range(len(metadata)):       
         transmission[ii][trans_wave>1215.67*(metadata[ii]['Z']+1)]=1.0
-              
+
     if(args.dla=='file'):
-        log.info('Adding DLAs from transmision file')
+        log.info('Adding DLAs from transmission file')
         min_trans_wave=np.min(trans_wave/1215.67 - 1)
         for ii in range(len(metadata)):
             if min_trans_wave < metadata[ii]['Z']: 
@@ -231,18 +231,19 @@ def simulate_one_healpix(ifilename,args,model,obsconditions,decam_and_wise_filte
 
     elif(args.dla=='random'):
         log.info('Adding DLAs randomly')
+        random_state_just_for_dlas = np.random.RandomState(seed)
         min_trans_wave=np.min(trans_wave/1215.67 - 1)
         for ii in range(len(metadata)):
             if min_trans_wave < metadata[ii]['Z']: 
-               idd=metadata['MOCKID'][ii]
-               dlass, dla_model = insert_dlas(trans_wave, metadata[ii]['Z'])
-               if len(dlass)>0:
-                  transmission[ii]=dla_model*transmission[ii]
-                  dla_z+=[idla['z'] for idla in dlass]
-                  dla_NHI+=[idla['N'] for idla in dlass]
-                  dla_id+=[idd]*len(dlass)    
+                idd=metadata['MOCKID'][ii]
+                dlass, dla_model = insert_dlas(trans_wave, metadata[ii]['Z'], rstate=random_state_just_for_dlas)
+                if len(dlass)>0:
+                    transmission[ii]=dla_model*transmission[ii]
+                    dla_z+=[idla['z'] for idla in dlass]
+                    dla_NHI+=[idla['N'] for idla in dlass]
+                    dla_id+=[idd]*len(dlass)
         log.info('Added {} DLAs'.format(len(dla_id)))
-
+        
     if args.dla is not None :
        dla_meta=Table()
        if len(dla_id)>0:    
@@ -297,7 +298,9 @@ def simulate_one_healpix(ifilename,args,model,obsconditions,decam_and_wise_filte
     if (args.balprob):
        if(args.balprob<=1. and args.balprob >0):
           log.info("Adding BALs with probability {}".format(args.balprob))
+          rnd_state = np.random.get_state() # save current random state
           tmp_qso_flux,meta_bal=bal.insert_bals(tmp_qso_wave,tmp_qso_flux, metadata['Z'], balprob=args.balprob,seed=seed)
+          np.random.set_state(rnd_state) # restore random state to get the same random numbers later as when we don't insert BALs
        else:
           log.error("Probability to add BALs is not between 0 and 1")
           sys.exit(1)
@@ -378,8 +381,7 @@ def simulate_one_healpix(ifilename,args,model,obsconditions,decam_and_wise_filte
         fibermap_columns={"MAG":mags}
     else :
         fibermap_columns=None
-
-    sim_spectra(qso_wave,qso_flux, args.program, obsconditions=obsconditions,spectra_filename=ofilename,sourcetype="qso", skyerr=args.skyerr,ra=metadata["RA"],dec=metadata["DEC"],targetid=targetid,meta=meta,seed=seed,fibermap_columns=fibermap_columns)
+    sim_spectra(qso_wave,qso_flux, args.program, obsconditions=obsconditions,spectra_filename=ofilename,sourcetype="qso", skyerr=args.skyerr,ra=metadata["RA"],dec=metadata["DEC"],targetid=targetid,meta=meta,seed=seed,fibermap_columns=fibermap_columns,use_poisson=False) # use Poisson = False to get reproducible results.
 
     if args.zbest is not None :
         log.info("Read fibermap")

py/desisim/scripts/quickspectra.py

@@ -24,7 +24,7 @@
 import matplotlib.pyplot as plt
 
 def sim_spectra(wave, flux, program, spectra_filename, obsconditions=None,
-                sourcetype=None, targetid=None, redshift=None, expid=0, seed=0, skyerr=0.0, ra=None, dec=None, meta=None, fibermap_columns=None):
+                sourcetype=None, targetid=None, redshift=None, expid=0, seed=0, skyerr=0.0, ra=None, dec=None, meta=None, fibermap_columns=None, fullsim=False,use_poisson=True):
     """
     Simulate spectra from an input set of wavelength and flux and writes a FITS file in the Spectra format that can
     be used as input to the redshift fitter.
@@ -49,7 +49,8 @@ def sim_spectra(wave, flux, program, spectra_filename, obsconditions=None,
         dec : numpy array with targets Dec (deg)
         meta : dictionnary, saved in primary fits header of the spectra file 
         fibermap_columns : add these columns to the fibermap
-    
+        fullsim : if True, write full simulation data in extra file per camera
+        use_poisson : if False, do not use numpy.random.poisson to simulate the Poisson noise. This is useful to get reproducible random realizations.
     """
     log = get_logger()
 
@@ -179,7 +180,7 @@ def sim_spectra(wave, flux, program, spectra_filename, obsconditions=None,
         psfconvolve=True)
 
     random_state = np.random.RandomState(seed)
-    sim.generate_random_noise(random_state)
+    sim.generate_random_noise(random_state,use_poisson=use_poisson)
 
     scale=1e17
     specdata = None
@@ -190,7 +191,14 @@ def sim_spectra(wave, flux, program, spectra_filename, obsconditions=None,
         resolution[camera.name] = np.tile(R.to_fits_array(), [nspec, 1, 1])
 
     skyscale = skyerr * random_state.normal(size=sim.num_fibers)
-
+
+    if fullsim :
+        for table in sim.camera_output :
+            band  = table.meta['name'].strip()[0]
+            table_filename=spectra_filename.replace(".fits","-fullsim-{}.fits".format(band))
+            table.write(table_filename,format="fits",overwrite=True)
+            print("wrote",table_filename)
+
     for table in sim.camera_output :
 
         wave = table['wavelength'].astype(float)
@@ -249,6 +257,7 @@ def parse(options=None):
     parser.add_argument('--seed', type=int, default=0, help="Random seed")
     parser.add_argument('--source-type', type=str, default=None, help="Source type (for fiber loss), among sky,elg,lrg,qso,bgs,star")
     parser.add_argument('--skyerr', type=float, default=0.0, help="Fractional sky subtraction error")
+    parser.add_argument('--fullsim',action='store_true',help="write full simulation data in extra file per camera, for debugging")
 
     if options is None:
         args = parser.parse_args()
@@ -345,5 +354,5 @@ def main(args=None):
 
     sim_spectra(input_wave, input_flux, args.program, obsconditions=obsconditions,
         spectra_filename=args.out_spectra,seed=args.seed,sourcetype=sourcetype,
-        skyerr=args.skyerr)
+        skyerr=args.skyerr,fullsim=args.fullsim)
 

py/desisim/spec_qa/high_level.py

@@ -84,7 +84,7 @@ def main():
     # Write
     summ_tab.meta = meta
     summ_tab.write(summ_file,format='ascii.ecsv',
-        formats=dict(MED_DZ='%8.6f',EFF='%5.3f',CAT_RATE='%6.4f'))#,clobber=True)
+        formats=dict(MED_DZ='%8.6f',EFF='%5.3f',CAT_RATE='%6.4f'))#,overwrite=True)
     '''
 
     # QA Figures

py/desisim/templates.py

@@ -222,14 +222,16 @@ def spectrum(self, oiiihbeta=None, oiihbeta=None, niihbeta=None,
         # Finally build the emission-line spectrum
         log10sigma = linesigma/LIGHT/np.log(10) # line-width [log-10 Angstrom]
         emspec = np.zeros(len(self.log10wave))
+
         for ii in range(len(line)):
             amp = line['flux'][ii] / line['wave'][ii] / np.log(10) # line-amplitude [erg/s/cm2/A]
             thislinewave = np.log10(line['wave'][ii] * (1.0+zshift))
             line['amp'][ii] = amp / (np.sqrt(2.0 * np.pi) * log10sigma)  # [erg/s/A]
 
             # Construct the spectrum [erg/s/cm2/A, rest]
-            emspec += amp * np.exp(-0.5 * (self.log10wave-thislinewave)**2 / log10sigma**2) \
-                      / (np.sqrt(2.0 * np.pi) * log10sigma)
+            jj = np.abs(self.log10wave-thislinewave) < 6*log10sigma
+            emspec[jj] += amp * np.exp(-0.5 * (self.log10wave[jj]-thislinewave)**2 / log10sigma**2) \
+                          / (np.sqrt(2.0 * np.pi) * log10sigma)
 
         return emspec, 10**self.log10wave, line
 

py/desisim/test/test_io.py

@@ -54,8 +54,10 @@ def test_simdir(self):
         self.assertTrue(x.endswith(str(expid)))
         x = io.simdir(int(night), expid)
         self.assertTrue(x.endswith(str(expid)))
-        x = io.simdir(night, mkdir=True)
-        self.assertTrue(os.path.exists(x))
+
+        #- If providing night, must also provide expid
+        with self.assertRaises(ValueError):
+            x = io.simdir(night, mkdir=True)
 
     def test_findfile(self):
         night = '20150102'