Add image weighting to mock observation #75
Add image weighting to mock observation #75nikos-triantafyllou wants to merge 2 commits into21cmfast:better_lc_noisefrom
Conversation
Added image weighting parameter to mock observation function and implemented weighting application for imaging.
nikos-triantafyllou
changed the title
for more information, see https://pre-commit.ci
| If an int is provided, all chunks will be | ||
| separated by the same number of slices. | ||
| image_weighting : str, optional | ||
| Filters the last image according to the chosen filter |
There was a problem hiding this comment.
This wording was slightly confusing to me. By "last" I think you mean the "ultimate" image, but "last" sounds like the last one in the array (like, the last realization). I would perhaps phrase it like "A weighting scheme for the transform from UV to image space, one of "blah", "blah", "blah"
| lightcone = intrinsic_lightcone_uv_nu | ||
| thermal_rms_uv = thermal_rms_uv_nu | ||
|
|
||
| if weight_type == "none" or weight_type == None: |
There was a problem hiding this comment.
nah, just take None. Don't be like matplotlib.
|
|
||
| elif weight_type in ( | ||
| "wiener", | ||
| "w", |
There was a problem hiding this comment.
only accept "wiener". Make people be explicit
| if weight_type == "none" or weight_type == None: | ||
| return lc_uv_nu | ||
|
|
||
| if weight_type in ("inverse_variance", "iv"): |
There was a problem hiding this comment.
only take "inverse_variance"
|
|
||
| if weight_type in ("inverse_variance", "iv"): | ||
| with np.errstate(divide="ignore", invalid="ignore"): | ||
| w = 1.0 / (2 * (thermal_rms_uv**2).value) # inverse variance weights |
There was a problem hiding this comment.
don't need the 2 here because it will get normalized out anyway.
| elif weight_type in ( | ||
| "wiener", | ||
| "w", | ||
| ): # essentially matching the power spectrum of the observed with the signal | ||
| # Estimate signal power spectrum from the measured signal and the known noise | ||
| # Wiener filter is defined usually for the power spectrum as w = (p_signal)/(p_signal + p_noise) | ||
| # the sigma in our case is applied in the real and imaginary parts separately. | ||
| # So the w computed should be further processed to w = np.sqrt(w)/2 (like in generating ICs for cosmological simulations) | ||
| p_signal = ( | ||
| np.abs(np.fft.rfft2(lightcone.value, axes=(0, 1))) ** 2 | ||
| ) # shape (nx, ny_rfft) | ||
|
|
||
| # Now p_signal is the pure signal, just fourier transformed. Let's take the mean along the frequency axis and repeat it to be used as a prior | ||
| p_signal = np.mean(p_signal, axis=2) | ||
|
|
||
| p_signal = p_signal[..., np.newaxis] | ||
| print(p_signal.shape) | ||
| p_signal = np.repeat(p_signal, 200, axis=2) | ||
| # Perhaps smooth it (maybe good but also sort of prior dependent) | ||
| # p_signal = gaussian_filter(p_signal, sigma=2) # from scipy.ndimage | ||
| p_noise = ( | ||
| 2 * (thermal_rms_uv.value) ** 2 | ||
| ) # our existing noise variance (real and imaginary parts so the variance is 2 sigma^2, hence the "2") | ||
| with np.errstate(divide="ignore", invalid="ignore"): | ||
| w = p_signal / (p_signal + p_noise) # Wiener filter, values in [0,1] | ||
| w = np.sqrt(w / 2) # applied on the 2D UV space on real and imag parts. | ||
|
|
There was a problem hiding this comment.
I think it is essentially always a bad idea to do this. In fact, doing essentially this is what caused the PAPER experiment (HERA precursor) to retract their first upper limit paper.
The problem is using the data itself to weight itself. This causes all kinds of weird statistical oddities. Instead, I would require the user to pass in a signal power spectrum they want to use to weight by. How they come up with that spectrum is up to them.
There was a problem hiding this comment.
Hmm, I see. I thought it was better to weigh from the data. Do you have an opinion on what would be the best for my use case? (i.e., fixed astro params, field-level, want to preserve as much structure and make as easy for the NN as possible)
There was a problem hiding this comment.
I.e., One cannot calculate the signal PS to use in the Wiener filter with cross-power spectra as you did with the noise PS calculation in the match_lc_noise_power_to_sense.ipynb tutorial?
There was a problem hiding this comment.
I am not sure I understand your connecting this to the cross-power of the noise I used in the matching to 21cmSense. Can you clarify what you mean there?
If it is a fixed astro params scenario, I would simply use the known isotropic power spectrum (averaged over many realizations, say) to weight the modes in your case.
2 participants
Added image weighting parameter to mock observation function and implemented weighting application for imaging.
Needs a bit more thinking before merging. The simple Wiener now is unrealistic, and the realistic Wiener is not good enough.
However, we should be able to do something similar with the cross power taking multiple realizations as @steven-murray did in the tutorials when comparing power spectra.
Here I show some demonstrations for the effects of the new function.
