Merge branch 'master' into specsim_streamline

.travis.yml

@@ -43,8 +43,8 @@ env:
         - ASTROPY_VERSION=1.1.1
         - SPHINX_VERSION=1.3
         - DESIUTIL_VERSION=1.3.3
-        - SPECTER_VERSION=0.3
-        - DESISPEC_VERSION=0.4
+        - SPECTER_VERSION=0.4
+        - DESISPEC_VERSION=0.4.1
         - DESITARGET_VERSION=0.3.2
         - DESIMODEL_VERSION=trunk
         - CONDA_INSTALL='conda install -c astropy-ci-extras --yes'

bin/pixsim-desi

@@ -46,6 +46,8 @@ parser.add_option("--trimxy", action="store_true", help="Trim image to fit spect
 parser.add_option("--randseed", type=int, help="random number seed")
 parser.add_option("--nspec", type=int, help="Number of spectra to simulate per camera [%default]", default=500)
 parser.add_option("--ncpu",  type=int, help="Number of cpu cores to use [%default]", default=multiprocessing.cpu_count() // 2)
+parser.add_option("--wavemin",  type=float, help="Minimum wavelength to simulate")
+parser.add_option("--wavemax",  type=float, help="Maximum wavelength to simulate")
 
 opts, args = parser.parse_args()
 

py/desisim/_version.py

@@ -1 +1 @@
-__version__ = '0.8.2.dev376'
+__version__ = '0.9.dev463'

py/desisim/pixsim.py

@@ -5,33 +5,33 @@
 import sys
 import os
 import os.path
-import time
-### import multiprocessing as mp
 
 import numpy as np
-import yaml
-from astropy.io import fits
 
 import desimodel.io
 import desispec.io
 from desispec.image import Image
+import desispec.cosmics
 
 from desisim import obs, io
 
 def simulate(night, expid, camera, nspec=None, verbose=False, ncpu=None,
-    trimxy=False, cosmics=None):
+    trimxy=False, cosmics=None, wavemin=None, wavemax=None):
     """
     Run pixel-level simulation of input spectra
     
     Args:
         night (string) : YEARMMDD
         expid (integer) : exposure id
         camera (str) : e.g. b0, r1, z9
+
+    Optional:
         nspec (int) : number of spectra to simulate
         verbose (boolean) : if True, print status messages
         ncpu (int) : number of CPU cores to use in parallel
         trimxy (boolean) : trim image to just pixels with input signal
         cosmics (str) : filename with dark images with cosmics to add
+        wavemin, wavemax (float) : min/max wavelength range to simulate
 
     Reads:
         $DESI_SPECTRO_SIM/$PIXPROD/{night}/simspec-{expid}.fits
@@ -77,6 +77,14 @@ def simulate(night, expid, camera, nspec=None, verbose=False, ncpu=None,
 
     phot = phot[ii]
 
+    #- Trim wavelenths if needed
+    if wavemin is not None:
+        ii = (wave >= wavemin)
+        phot = phot[:, ii]
+    if wavemax is not None:
+        ii = (wave <= wavemax)
+        phot = phot[:, ii]
+
     #- check if simulation has less than 500 input spectra
     if phot.shape[0] < nspec:
         nspec = phot.shape[0]
@@ -114,22 +122,22 @@ def simulate(night, expid, camera, nspec=None, verbose=False, ncpu=None,
         cosmics = io.read_cosmics(cosmics, expid, shape=img.shape)
         pix = np.random.poisson(img) + cosmics.pix
         readnoise = cosmics.meta['RDNOISE']
-        ivar = 1.0/(pix.clip(0) + readnoise**2)
-        mask = cosmics.mask
     #- Or just add noise
     else:
-        params = desimodel.io.load_desiparams()
         channel = camera[0].lower()
         readnoise = params['ccd'][channel]['readnoise']
         pix = np.random.poisson(img) + np.random.normal(scale=readnoise, size=img.shape)
-        ivar = 1.0/(pix.clip(0) + readnoise**2)
-        mask = np.zeros(img.shape, dtype=np.uint16)
-
+
+    ivar = 1.0/(pix.clip(0) + readnoise**2)
+    mask = np.zeros(img.shape, dtype=np.uint16)
+
     #- Metadata to be included in pix file header is in the fibermap header
     #- TODO: this is fragile; consider updating fibermap to use astropy Table
     #- that includes the header rather than directly assuming FITS as the
     #- underlying format.
+    simdir = os.path.dirname(simfile)
     fibermapfile = desispec.io.findfile('fibermap', night=night, expid=expid)
+    fibermapfile = os.path.join(simdir, os.path.basename(fibermapfile))
     fm, fmhdr = desispec.io.read_fibermap(fibermapfile, header=True)
     meta = dict()
     try:
@@ -146,7 +154,13 @@ def simulate(night, expid, camera, nspec=None, verbose=False, ncpu=None,
     meta['AIRMASS'] = simspec.header['AIRMASS']
 
     image = Image(pix, ivar, mask, readnoise=readnoise, camera=camera, meta=meta)
+
+    #- In-place update of the image cosmic ray mask
+    if cosmics is not None:
+        desispec.cosmics.reject_cosmic_rays(image)
+
     pixfile = desispec.io.findfile('pix', night=night, camera=camera, expid=expid)
+    pixfile = os.path.join(simdir, os.path.basename(pixfile))
     desispec.io.write_image(pixfile, image)
 
     if verbose:

py/desisim/targets.py

@@ -2,6 +2,8 @@
 Utility functions for working with simulated targets
 """
 
+from __future__ import absolute_import, division, print_function
+
 import os
 import sys
 import string
@@ -14,6 +16,7 @@
 
 from desispec import brick
 from desispec.io.fibermap import empty_fibermap
+from desitarget.targetmask import desi_mask, bgs_mask, mws_mask
 
 from desisim import io
 
@@ -36,6 +39,7 @@ def sample_objtype(nobj, flavor):
       LRGs and QSOs in early passes and more ELGs in later passes.
     - Also ensures at least 2 sky and 1 stdstar, even if nobj is small
     """
+    flavor = flavor.upper()
 
     #- Load target densities
     #- TODO: what about nobs_boss (BOSS-like LRGs)?
@@ -44,12 +48,6 @@ def sample_objtype(nobj, flavor):
     tgt = yaml.load(fx)
     fx.close()
 
-    # DARK or BRIGHT have a combination of targets
-    if string.lower(flavor) == 'dark':
-        ntgt = float(tgt['nobs_lrg'] + tgt['nobs_elg'] + tgt['nobs_qso'] + tgt['nobs_lya'] + tgt['ntarget_badqso'])
-    elif string.lower(flavor) == 'bright':
-        ntgt = float(tgt['nobs_BG'] + tgt['nobs_MWS'])
-
     # initialize so we can ask for 0 of some kinds of survey targets later
     nlrg = nqso = nelg = nmws = nbgs = nbgs = nmws = 0
     #- Fraction of sky and standard star targets is guaranteed
@@ -67,39 +65,46 @@ def sample_objtype(nobj, flavor):
     nsci = nobj - (nsky+nstd)
     true_objtype  = ['SKY']*nsky + ['STD']*nstd
 
-    if (string.lower(flavor) == 'mws'):
+    if (flavor == 'MWS'):
         true_objtype  +=  ['MWS_STAR']*nsci
-    elif (string.lower(flavor) == 'qso'):
+    elif (flavor == 'QSO'):
         true_objtype  +=  ['QSO']*nsci
-    elif (string.lower(flavor) == 'elg'):
+    elif (flavor == 'ELG'):
         true_objtype  +=  ['ELG']*nsci    
-    elif (string.lower(flavor) == 'lrg'):
+    elif (flavor == 'LRG'):
         true_objtype  +=  ['LRG']*nsci
-    elif (string.lower(flavor) == 'std'):
+    elif (flavor == 'STD'):
         true_objtype  +=  ['STD']*nsci
-    elif (string.lower(flavor) == 'bgs'):
+    elif (flavor == 'BGS'):
         true_objtype  +=  ['BGS']*nsci
-    elif (string.lower(flavor) == 'bright'):
+    elif (flavor == 'BRIGHT'):
         #- BGS galaxies and MWS stars
-        nbgs = min(nsci, np.random.poisson(nsci * tgt['nobs_BG'] / ntgt))
-        nmws = nsci - nbgs
+        ntgt = float(tgt['nobs_BG'] + tgt['nobs_MWS'])
+        prob_bgs = tgt['nobs_BG'] / ntgt
+        prob_mws = 1 - prob_bgs
+
+        p = [prob_bgs, prob_mws]
+        nbgs, nmws = np.random.multinomial(nsci, p)
+
         true_objtype += ['BGS']*nbgs
         true_objtype += ['MWS_STAR']*nmws
-    elif (string.lower(flavor) == 'dark'):
+    elif (flavor == 'DARK'):
         #- LRGs ELGs QSOs
-        nlrg = np.random.poisson(nsci * tgt['nobs_lrg'] / ntgt)
-
-        nqso = np.random.poisson(nsci * (tgt['nobs_qso'] + tgt['nobs_lya']) / ntgt)
-        nqso_bad = np.random.poisson(nsci * (tgt['ntarget_badqso']) / ntgt)
-
-        nelg = nobj - (nlrg+nqso+nqso_bad+nsky+nstd)
+        ntgt = float(tgt['nobs_lrg'] + tgt['nobs_elg'] + tgt['nobs_qso'] + tgt['nobs_lya'] + tgt['ntarget_badqso'])
+        prob_lrg = tgt['nobs_lrg'] / ntgt
+        prob_elg = tgt['nobs_elg'] / ntgt
+        prob_qso = (tgt['nobs_qso'] + tgt['nobs_lya']) / ntgt
+        prob_badqso = tgt['ntarget_badqso'] / ntgt
+
+        p = [prob_lrg, prob_elg, prob_qso, prob_badqso]
+        nlrg, nelg, nqso, nbadqso = np.random.multinomial(nsci, p)
 
         true_objtype += ['ELG']*nelg
         true_objtype += ['LRG']*nlrg
-        true_objtype += ['QSO']*nqso + ['QSO_BAD']*nqso_bad
+        true_objtype += ['QSO']*nqso + ['QSO_BAD']*nbadqso
     else:
         raise ValueError("Do not know the objtype mix for flavor "+ flavor)
-
+        
     assert len(true_objtype) == nobj, \
         'len(true_objtype) mismatch for flavor {} : {} != {}'.format(\
         flavor, len(true_objtype), nobj)
@@ -132,6 +137,8 @@ def get_targets(nspec, flavor, tileid=None):
     else:
         tile_ra, tile_dec = io.get_tile_radec(tileid)
 
+    flavor = flavor.upper()
+
     #- Get distribution of target types
     true_objtype, target_objtype = sample_objtype(nspec, flavor)
 
@@ -167,27 +174,33 @@ def get_targets(nspec, flavor, tileid=None):
 
         # Simulate spectra
         if objtype == 'SKY':
+            fibermap['DESI_TARGET'][ii] = desi_mask.SKY
             continue
 
         elif objtype == 'ELG':
             from desisim.templates import ELG
             elg = ELG(wave=wave)
             simflux, wave1, meta = elg.make_templates(nmodel=nobj)
+            fibermap['DESI_TARGET'][ii] = desi_mask.ELG
 
         elif objtype == 'LRG':
             from desisim.templates import LRG
             lrg = LRG(wave=wave)
             simflux, wave1, meta = lrg.make_templates(nmodel=nobj)
+            fibermap['DESI_TARGET'][ii] = desi_mask.LRG
 
         elif objtype == 'BGS':
             from desisim.templates import BGS
             bgs = BGS(wave=wave)
             simflux, wave1, meta = bgs.make_templates(nmodel=nobj)
+            fibermap['DESI_TARGET'][ii] = desi_mask.BGS_ANY
+            fibermap['BGS_TARGET'][ii] = bgs_mask.BGS_BRIGHT
 
         elif objtype == 'QSO':
             from desisim.templates import QSO
             qso = QSO(wave=wave)
             simflux, wave1, meta = qso.make_templates(nmodel=nobj)
+            fibermap['DESI_TARGET'][ii] = desi_mask.QSO
 
         # For a "bad" QSO simulate a normal star without color cuts, which isn't
         # right. We need to apply the QSO color-cuts to the normal stars to pull
@@ -196,6 +209,7 @@ def get_targets(nspec, flavor, tileid=None):
             from desisim.templates import STAR
             star = STAR(wave=wave)
             simflux, wave1, meta = star.make_templates(nmodel=nobj)
+            fibermap['DESI_TARGET'][ii] = desi_mask.QSO
 
         elif objtype == 'STD':
             from desisim.templates import STAR
@@ -204,12 +218,18 @@ def get_targets(nspec, flavor, tileid=None):
             gg = (16.0, 19.5)
             simflux, wave1, meta = \
                 star.make_templates(nmodel=nobj, rmagrange=rr, gmagrange=gg)
+            fibermap['DESI_TARGET'][ii] = desi_mask.STD_FSTAR
 
         elif objtype == 'MWS_STAR':
             from desisim.templates import STAR
             star = STAR(wave=wave)
             # todo: mag ranges for different flavors of STAR targets should be in desimodel
             simflux, wave1, meta = star.make_templates(nmodel=nobj,rmagrange=(15.0,20.0))
+            fibermap['DESI_TARGET'][ii] = desi_mask.MWS_ANY
+            fibermap['MWS_TARGET'][ii] = mws_mask.MWS_PLX  #- ???
+
+        else:
+            raise ValueError('Unable to simulate OBJTYPE={}'.format(objtype))
 
         truth['FLUX'][ii] = 1e17 * simflux
         truth['UNITS'] = '1e-17 erg/s/cm2/A'
@@ -232,7 +252,7 @@ def get_targets(nspec, flavor, tileid=None):
             truth['D4000'][ii] = meta['D4000']
             truth['VDISP'][ii] = meta['VDISP']
 
-        if objtype == 'BGS':
+        if objtype == 'BGS':            
             truth['HBETAFLUX'][ii] = meta['HBETAFLUX']
             truth['D4000'][ii] = meta['D4000']
             truth['VDISP'][ii] = meta['VDISP']
@@ -256,7 +276,6 @@ def get_targets(nspec, flavor, tileid=None):
     fibermap['TARGETID'] = np.random.randint(sys.maxint, size=nspec)
     fibermap['TARGETCAT'] = np.zeros(nspec, dtype='|S20')
     fibermap['LAMBDAREF'] = np.ones(nspec, dtype=np.float32)*5400
-    fibermap['TARGET_MASK0'] = np.zeros(nspec, dtype='i8')
     fibermap['RA_TARGET'] = ra
     fibermap['DEC_TARGET'] = dec
     fibermap['X_TARGET'] = x
@@ -433,7 +452,7 @@ def load_results(self, targets_file):
             self.assigned_type =  fin[1].data['ASSIGNEDTYPE']
         except Exception, e:
             import traceback
-            print 'ERROR in get_tiles'
+            print('ERROR in get_tiles')
             traceback.print_exc()
             raise e
 

py/desisim/test/test_targets.py

@@ -19,6 +19,13 @@ def test_get_targets(self):
         n = 5
         for flavor in ['DARK', 'BRIGHT', 'LRG', 'ELG', 'QSO', 'BGS', 'MWS']:
             fibermap, truth = desisim.targets.get_targets(n, flavor)
+            fibermap, truth = desisim.targets.get_targets(n, flavor.lower())
+
+    def test_sample_objtype(self):
+        for flavor in ['DARK', 'BRIGHT', 'LRG', 'ELG', 'QSO', 'BGS', 'MWS']:
+            for n in [5,10,50]:
+                for i in range(10):
+                    truetype, targettype = desisim.targets.sample_objtype(n, flavor)
 
 #- This runs all test* functions in any TestCase class in this file
 if __name__ == '__main__':