Smoothing source contribution to noise in quickquasars DESI mocks #566
Conversation
…oise for DESI mocks. The value is the Gaussian sigma in A. This is implemented with FFTs on a padded array.
py/desisim/scripts/quickspectra.py
Outdated
| # Smoothing source electron numbers only works for DESI mocks | ||
| if specsim_config_file != "eboss" and source_contribution_smoothing > 0: | ||
| _smooth_source_variance(sim.camera_output, source_contribution_smoothing, dwave_out) | ||
|
|
||
| random_state = np.random.RandomState(seed) | ||
| sim.generate_random_noise(random_state,use_poisson=use_poisson) |
There was a problem hiding this comment.
@julienguy , @p-slash , @alxogm - shouldn't we hack the estimate of the variance after we have generated the random noise?
There was a problem hiding this comment.
Yes. We should change the variance model after generating the noise.
| kernel_k = np.exp(-(kvals*sigmapix)**2/2.) | ||
| snumsource_k = np.fft.rfft(padded_arr, axis=0)*kernel_k[:, None] | ||
|
|
||
| return np.fft.irfft(snumsource_k, n=arrsize, axis=0)[pad_size:-pad_size] |
There was a problem hiding this comment.
I haven't checked the pythonology of arrays and FFT conventions here, but I guess you have checked this.
There was a problem hiding this comment.
It doesn't look like the kernel is properly normalized here.
As for the fft calls and padding, I am personally using scipy.signal.fftconvolve (see for instance https://github.com/desihub/desispec/blob/d925b926e9e3767d2cbd0457e3d8afa42b55ecfb/py/desispec/image_model.py#L132 )
There was a problem hiding this comment.
I can switch to fftconvolve. Any tricks for edge effects?
Why is this implementation not normalized though? kernel_k[0]=1, so sum of kernel_x=1. That's what I thought at least.
There was a problem hiding this comment.
(my bad, I was thinking about a kernel in configuration space, as it is the case in the argument of fftconvolve)
There was a problem hiding this comment.
A part from my comment on when to smooth the noise variance, I think it looks good. I would like to see a comparison of reported values of IVAR from different runs of the same spectrum, with different values of smoothing, just to check that indeed the reported values are smoothed out.
As I mentioned, I think the flux values should be exactly the same, including the noise realisation, and that only the reported values of noise variance should have changed.
|
Sorry, this is the first time I'm seeing this. This is a subtle point. If we were to have the exact flux and noise, but different reported variance, that would be wrong. That is just the wrong noise estimate for that realization. The point is to decouple the generated noise from the source flux, so they will be different for different smoothing scales. In this implementation, we smooth the variance first (only the source contribution) and generate noise from this smoothed variance. This way, the generated noise is decoupled from the source flux, and we still know the exact noise that is put in. Furthermore, using these reported variances as inverse weights will be unbiased (less biased) in statistics estimations. In terms of P1D and Xi3D estimation, I think we should still smooth the noise in data and mocks, but only in its contribution to weights. Noise power subtraction should use the exact reported variance to correctly remove the noise power. The difference is that here we are smoothing out all contributions to noise, whereas the implementation in quickquasars smoothes out only the source contribution. |
|
@julienguy , could you comment on this? I thought that what I suggested above is what the real pipeline does. |
|
@julienguy , @p-slash - we forgot to discuss this at the telecon yesterday. Should we have a quick call to make sure we all agree on this? Do you prefer to discuss this online? |
|
Also, @p-slash - can you move the branch to the desisim repo, instead of your own fork? César Ramírez would look at it in the context of BAO analyses. |
|
Not sure how. I can try. Can't he pull from my branch to his local repo? I would be fine with a chat or over text on Slack. I don't check here as often:) |
|
@p-slash , @julienguy - any update on this discussion? |
|
I commented we should change the variance model after generating the noise. Naim @p-slash I can make the change if you are not sure how to do this. |
|
I moved the function after sim.generate_random_noise(random_state,use_poisson=use_poisson). That should do the trick. I HAVE NOT tested it though. Sorry I've been busy with other things. |
|
@p-slash , @julienguy , @alxogm - what is the status of this PR? |
|
Sorry, I missed to follow this PR. Also I think we haven't tested the effect on the analysis. I can do early next week. |
|
Thanks Alma |
|
Hi @p-slash I have only requested to add one line to run log to inform the smoothing was applied. |
|
Given our discussion at the workshop, should we change the default for |
| @@ -189,6 +229,10 @@ def sim_spectra(wave, flux, program, spectra_filename, obsconditions=None, | |||
| random_state = np.random.RandomState(seed) | |||
| sim.generate_random_noise(random_state,use_poisson=use_poisson) | |||
|
|
|||
| # Smoothing source electron numbers only works for DESI mocks | |||
| if specsim_config_file != "eboss" and source_contribution_smoothing > 0: | |||
There was a problem hiding this comment.
Could we add a line to the run log to inform what smoothing was applied?
There was a problem hiding this comment.
@p-slash I suggest here you can add a line for the log
There is now an option in quickquasars to smooth out the source contribution to noise for DESI mocks in order to de-bias the weights in Lya forest analyses. The value passed will the Gaussian sigma for the smoothing in Angstrom. This smoothing is implemented with FFTs on a padded array in quickspectra.py.