No poisson

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,17 +231,19 @@ def simulate_one_healpix(ifilename,args,model,obsconditions,decam_and_wise_filte
 
     elif(args.dla=='random'):
         log.info('Adding DLAs randomly')
+        saved_rnd_state = np.random.get_state()
         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'], seed=np.random.randint(2**32)) # need to pass a random seed otherwise insert_dla breaks the random sequence
+                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)))
+        np.random.set_state(saved_rnd_state) # restore random state to get the same random numbers later as when we don't insert DLAs
 
     if args.dla is not None :
        dla_meta=Table()
@@ -297,7 +299,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 +382,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,8 @@ 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()
 
@@ -178,10 +179,10 @@ def sim_spectra(wave, flux, program, spectra_filename, obsconditions=None,
     sim = desisim.simexp.simulate_spectra(wave, flux, fibermap=frame_fibermap,
         obsconditions=obsconditions, redshift=redshift, seed=seed,
         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