Merge pull request #393 from desihub/no-poisson

No poisson

.travis.yml

@@ -53,7 +53,7 @@ env:
         - 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

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/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, fullsim=False):
+                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.
@@ -50,7 +50,7 @@ def sim_spectra(wave, flux, program, spectra_filename, obsconditions=None,
         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()
 
@@ -180,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

py/desisim/test/test_quickgen.py

@@ -1,9 +1,10 @@
 """
 Run unit tests through quickgen
 """
-import unittest, os
+import unittest, os, sys
 from uuid import uuid1
 from shutil import rmtree
+import subprocess
 
 from astropy.io import fits
 import numpy as np
@@ -84,6 +85,25 @@ def setUpClass(cls):
         else:
             cls.cosmics = None
 
+        #- Try to locate where the quickgen script will be
+        cls.topDir = os.path.dirname( # top-level
+            os.path.dirname( # build/
+                os.path.dirname( # lib/
+                    os.path.dirname( # desispec/
+                        os.path.dirname(os.path.abspath(__file__)) # test/
+                        )
+                    )
+                )
+            )
+        cls.binDir = os.path.join(cls.topDir,'bin')
+
+        if not os.path.isdir(cls.binDir):
+            cls.topDir = os.getcwd()
+            cls.binDir = os.path.join(cls.topDir, 'bin')
+
+        if not os.path.isdir(cls.binDir):
+            raise RuntimeError('Unable to auto-locate desisim/bin from {}'.format(__file__))
+
     def setUp(self):
         self.night = '20150105'
         self.expid = 124
@@ -227,13 +247,36 @@ def test_quickgen_seed_simspec(self):
         simspec2 = io.findfile('simspec', night, expid2)
         fibermap2 = desispec.io.findfile('fibermap', night, expid2)
 
+        s0=fits.open(simspec0)
+        s1=fits.open(simspec1)
+        s2=fits.open(simspec2)
+
+        self.assertEqual(s0['TRUTH'].data['OBJTYPE'][0], s1['TRUTH'].data['OBJTYPE'][0])
+        self.assertEqual(s0['TRUTH'].data['REDSHIFT'][0], s1['TRUTH'].data['REDSHIFT'][0])
+        self.assertNotEqual(s0['TRUTH'].data['REDSHIFT'][0], s2['TRUTH'].data['REDSHIFT'][0])
+
+        self.assertTrue(np.all(s0['WAVE'].data == s1['WAVE'].data))
+        self.assertTrue(np.all(s0['FLUX'].data == s1['FLUX'].data))
+
         # generate quickgen output for each exposure
-        cmd = "quickgen --simspec {} --fibermap {} --seed 1".format(simspec0,fibermap0)
-        quickgen.main(quickgen.parse(cmd.split()[1:]))
-        cmd = "quickgen --simspec {} --fibermap {} --seed 1".format(simspec1,fibermap1)
-        quickgen.main(quickgen.parse(cmd.split()[1:]))
-        cmd = "quickgen --simspec {} --fibermap {} --seed 2".format(simspec2,fibermap2)
-        quickgen.main(quickgen.parse(cmd.split()[1:]))
+        # spawn with subprocess so that it will really be restarting fresh;
+        # using quickgen.main re-uses a specsim Simulator that consumes random
+        # state when using it.
+        # See https://github.com/desihub/specsim/issues/94
+        cmd = "{} {}/quickgen --simspec {} --fibermap {} --seed 1".format(
+            sys.executable, self.binDir, simspec0, fibermap0)
+        ### quickgen.main(quickgen.parse(cmd.split()[1:]))
+        subprocess.check_call(cmd.split())
+
+        cmd = "{} {}/quickgen --simspec {} --fibermap {} --seed 1".format(
+            sys.executable, self.binDir, simspec1, fibermap1)
+        ### quickgen.main(quickgen.parse(cmd.split()[1:]))
+        subprocess.check_call(cmd.split())
+
+        cmd = "{} {}/quickgen --simspec {} --fibermap {} --seed 2".format(
+            sys.executable, self.binDir, simspec2, fibermap2)
+        ### quickgen.main(quickgen.parse(cmd.split()[1:]))
+        subprocess.check_call(cmd.split())
 
         cframe0=desispec.io.findfile("cframe",night,expid0,camera)
         cframe1=desispec.io.findfile("cframe",night,expid1,camera)
@@ -246,14 +289,6 @@ def test_quickgen_seed_simspec(self):
         self.assertTrue(np.all(cf0.ivar == cf1.ivar))
         self.assertTrue(np.any(cf0.flux != cf2.flux))   #- different seed
 
-        s0=fits.open(simspec0)
-        s1=fits.open(simspec1)
-        s2=fits.open(simspec2)
-
-        self.assertEqual(s0['TRUTH'].data['OBJTYPE'][0], s1['TRUTH'].data['OBJTYPE'][0])
-        self.assertEqual(s0['TRUTH'].data['REDSHIFT'][0], s1['TRUTH'].data['REDSHIFT'][0])
-        self.assertNotEqual(s0['TRUTH'].data['REDSHIFT'][0], s2['TRUTH'].data['REDSHIFT'][0])
-
     #- Test that higher airmass makes noisier spectra using simspec as input
     ### @unittest.skipIf('TRAVIS_JOB_ID' in os.environ, 'Skipping memory hungry quickgen/specsim test on Travis')
     def test_quickgen_airmass_simspec(self):
@@ -400,3 +435,10 @@ def test_quickgen_moonzenith_simspec(self):
 #- This runs all test* functions in any TestCase class in this file
 if __name__ == '__main__':
     unittest.main()
+
+def test_suite():
+    """Allows testing of only this module with the command::
+
+        python setup.py test -m desisim.test.test_io
+    """
+    return unittest.defaultTestLoader.loadTestsFromName(__name__)

py/desisim/test/test_quickspectra.py

@@ -55,10 +55,10 @@ def _check_spectra_match(self, sp1, sp2, invert=False):
             self.assertTrue(np.allclose(sp1.wave[x], sp2.wave[x]))
             if invert:
                 self.assertFalse(np.allclose(sp1.flux[x], sp2.flux[x]))
-                self.assertFalse(np.allclose(sp1.ivar[x], sp2.ivar[x]))
             else:
                 self.assertTrue(np.allclose(sp1.flux[x], sp2.flux[x]))
                 self.assertTrue(np.allclose(sp1.ivar[x], sp2.ivar[x]))
+
 
     def test_quickspectra(self):
         cmd = 'quickspectra -i {} -o {} --seed 1'.format(self.inspec_txt, self.outspec1)