diff --git a/.travis.yml b/.travis.yml index 6bd2f1507..8b70cc5cc 100644 --- a/.travis.yml +++ b/.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 diff --git a/py/desisim/dla.py b/py/desisim/dla.py index aaea23ba6..06366a821 100644 --- a/py/desisim/dla.py +++ b/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) diff --git a/py/desisim/scripts/quickquasars.py b/py/desisim/scripts/quickquasars.py index 30db505d5..af65e8308 100644 --- a/py/desisim/scripts/quickquasars.py +++ b/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") diff --git a/py/desisim/scripts/quickspectra.py b/py/desisim/scripts/quickspectra.py index 65d06c7af..575b78341 100644 --- a/py/desisim/scripts/quickspectra.py +++ b/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 diff --git a/py/desisim/test/test_quickgen.py b/py/desisim/test/test_quickgen.py index e26b0c6d8..5eca8e90d 100644 --- a/py/desisim/test/test_quickgen.py +++ b/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__) diff --git a/py/desisim/test/test_quickspectra.py b/py/desisim/test/test_quickspectra.py index 78a7e9747..7c135ff50 100644 --- a/py/desisim/test/test_quickspectra.py +++ b/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)