No poisson #393

julienguy · 2018-07-20T23:31:04Z

numpy.random.poisson has the drawback of using a varying number of random numbers depending on the mean value (the argument of the function). this is problematic for some applications like testing the effect of DLA on the Lya correlation function where one would like to preserve the random numbers used to simulate the instrumental noise when adding DLAs or not. This PR fixes this by allowing the use of numpy.random.normal instead of numpy.random.poisson in specsim.simulator.Simulator.generate_random_noise.

Solves issue #386.
(Need to rerun tests when the equivalent PR in specsim is merged.)

…dom seauence, use random.normal insteal of random.poisson to simulate the noise to preserve random numbers per wavelength bin

dkirkby · 2018-07-21T08:10:30Z

py/desisim/scripts/quickspectra.py

@@ -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.


-> "This is useful to decouple the (normalized) noise contributions from the predicted mean signal."

dkirkby · 2018-07-21T08:11:28Z

py/desisim/scripts/quickspectra.py

@@ -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)
-
+    


invisible white space (https://www.python.org/dev/peps/pep-0008/?#id28)

dkirkby · 2018-07-21T08:13:28Z

py/desisim/scripts/quickquasars.py

@@ -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()


Its much easier to manage different parts of a program consuming random numbers independently if they each have their own private np.random.RandomState. A corollary is that you should avoid using the global RandomState with direct calls to np.random....

dkirkby · 2018-07-21T08:14:53Z

py/desisim/scripts/quickquasars.py

-                  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


Pass in a dla_random_state here that was created during initialization (instead of a seed).

dkirkby · 2018-07-21T08:15:17Z

py/desisim/scripts/quickquasars.py

        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


This is not necessary if the DLA code has its own RandomState object.

yes, now that I understand where the change of noise is coming from.

dkirkby · 2018-07-21T08:15:37Z

py/desisim/scripts/quickquasars.py

@@ -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


(Similar comments here)

…no need to restore global rnd state 'index' after dla insertion

dkirkby · 2018-07-24T07:39:16Z

I am restarting the travis tests after merging desihub/specsim#92. If those pass, this looks ready to merge.

dkirkby · 2018-07-24T07:50:22Z

It looks like travis is still stuck on a docstring formatting error:

Warning, treated as error:
/home/travis/build/desihub/desisim/py/desisim/scripts/quickspectra.py:
  docstring of desisim.scripts.quickspectra.sim_spectra:25:Unexpected indentation.

dkirkby · 2018-07-24T08:01:53Z

Taking a closer look at this package's travis config, tests are hardcoded to use specsim v0.11 so are not using the update to specsim master. I will tag a new version of specsim soon (with desimodel updates) so I suggest we wait for that before updating the travis config here, but this PR could be merged now if its blocking someone.

andreufont · 2018-07-24T08:12:31Z

@dkirkby - It is not critical, but it would be great to have this solved soonish, since it will allow several tests in the Lyman alpha working group.

dkirkby · 2018-07-26T15:51:04Z

I just tagged specsim v0.12 and updated this package's travis tests to use this new version. If they pass this is ready to merge & tag, then both tags can be included in the 18.7 release.

andreufont · 2018-07-26T15:54:56Z

That's great! Fingers crossed :-)

dkirkby · 2018-07-26T16:27:04Z

Travis tests are still failing, but no longer due to specsim (as far as I can tell).

@andreufont You should be able to use specsim master with this desisim branch for the Lya studies now, if you don't want to wait until this propagates into the 18.7 release.

julienguy · 2018-07-26T17:22:08Z

desisim docstring issue ... I am having a look at this

julienguy · 2018-07-26T17:50:00Z

a curious test failure that should never have succeeded: in quickspectra, when the seed is changed the flux in spectra is expected to change but not the ivar. however the test is requiring a different ivar (??). I am changing this.

julienguy · 2018-07-26T19:36:35Z

quickgen test failing (nothing to do with quickspectra or quickquasars ...)

sbailey · 2018-07-26T20:29:22Z

I confirm that the quickgen tests passes with specsim v0.11 but fails with specsim v0.12. Now to dig deeper...

julienguy · 2018-07-26T21:02:24Z

test don't pass because quickgen doesn't give the same result when called twice with the same random seed.
specsim.instrument.initialize(config) has a np.random.RandomState() if config.instrument.offset has no attribute 'seed'. this is what breaks the random number sequence. Setting a seed value in a config file wouldn't make sense, so one would need to alter the config object before it is used in the constructor of instrument ... I am not sure where one should do that and I don't think we should have a random seed in a configuration object anyways.
a short term fix would be load the config, modify it, then create the simulator object in desisim.specsim.get_simulator(with a new argument seed) ... but that looks like a band-aid solution.

…andom seeds

sbailey · 2018-07-26T21:55:39Z

I spun the specsim part to desihub/specsim#94 . This quickgen test was using a seed for generating the instrument offsets random state, but since it was using the same Simulator object multiple times in a row, its results depended upon this history of what else was simulated first. i.e. simulating the same spectra twice in a row using the same simulator doesn't give the same answer.

For now I've "fixed" this in the desisim tests by spawning quickgen using subprocess.call instead of quickgen.main() so that I won't re-use the same specsim.Simulator on the second call.

sbailey · 2018-07-26T22:11:45Z

Arg. My subprocess.call fix works when desisim/bin is in $PATH, but not for the bare-bones testing environment pre-installation and without setting $PATH. Working on it...

julienguy added 3 commits July 20, 2018 15:55

add a comment about seeds

3adb616

add option use_poisson (with default=True)

7bc9c30

save and restore random states, use seed for DLA to keep coherent ran…

2ce432d

…dom seauence, use random.normal insteal of random.poisson to simulate the noise to preserve random numbers per wavelength bin

julienguy requested a review from dkirkby July 20, 2018 23:47

dkirkby requested changes Jul 21, 2018

View reviewed changes

julienguy added 2 commits July 23, 2018 09:15

two white spaces

50a8896

create the rnd state 'instance' here rather than in the dla routine. …

3968d41

…no need to restore global rnd state 'index' after dla insertion

dkirkby approved these changes Jul 23, 2018

View reviewed changes

Use specsim v0.12 for travis tests

31d4fcf

docstring indentation

07e8f5e

modify ivar test

ebb612f

sbailey mentioned this pull request Jul 26, 2018

specsim.instrument and config random state desihub/specsim#94

Open

spawn quickgen with subprocess instead of quickgen.main for testing r…

1d96d09

…andom seeds

find desisim/bin path

0d7ca01

sbailey merged commit b119762 into master Jul 26, 2018

sbailey deleted the no-poisson branch July 26, 2018 23:09

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

No poisson #393

No poisson #393

julienguy commented Jul 20, 2018

dkirkby Jul 21, 2018

dkirkby Jul 21, 2018

julienguy Jul 23, 2018

dkirkby Jul 21, 2018

dkirkby Jul 21, 2018

julienguy Jul 23, 2018

dkirkby Jul 21, 2018

julienguy Jul 23, 2018

dkirkby Jul 21, 2018

dkirkby commented Jul 24, 2018

dkirkby commented Jul 24, 2018

dkirkby commented Jul 24, 2018

andreufont commented Jul 24, 2018

dkirkby commented Jul 26, 2018

andreufont commented Jul 26, 2018

dkirkby commented Jul 26, 2018

julienguy commented Jul 26, 2018

julienguy commented Jul 26, 2018

julienguy commented Jul 26, 2018

sbailey commented Jul 26, 2018

julienguy commented Jul 26, 2018

sbailey commented Jul 26, 2018

sbailey commented Jul 26, 2018

		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

No poisson #393

No poisson #393

Conversation

julienguy commented Jul 20, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dkirkby commented Jul 24, 2018

dkirkby commented Jul 24, 2018

dkirkby commented Jul 24, 2018

andreufont commented Jul 24, 2018

dkirkby commented Jul 26, 2018

andreufont commented Jul 26, 2018

dkirkby commented Jul 26, 2018

julienguy commented Jul 26, 2018

julienguy commented Jul 26, 2018

julienguy commented Jul 26, 2018

sbailey commented Jul 26, 2018

julienguy commented Jul 26, 2018

sbailey commented Jul 26, 2018

sbailey commented Jul 26, 2018