Pk1d diff noise estimators #997

corentinravoux · 2023-04-26T10:25:23Z

Short pull request implementing the latest tests I realized on the diff estimator for EDR. We forgot to report them from the divergent Pk1d_DESI_changes branch.
I tested different methods and I put the one which is now used for the EDR paper as default.

iprafols

Looks fine to me. Most of the changes are just formatting changes. If I understood it correctly, the only algorithmic changes are those of exp_dif_desi. Maybe someone more expert on Pk1d can review it to make sure it's ok

Waelthus

The addition of the different methods used is great!
It would've been nice to start on an up-to-date version of the file when implementing this, now we'd either need to go back to old variable names or need to fix the variable names in the PR to match conventions.
There's some suggestions regarding re-ordering, potentially removing some of the for-loops to shorten the code overhead, and some questions regarding the method used further up.

Waelthus · 2023-04-27T08:43:33Z

py/picca/pk1d/prep_pk1d.py

+    mask_targetid,
+    method_alpha="desi_array",
+    use_only_even=False,
+):
    """Compute the difference between exposures.


Please update the docstring to the new arguments. If file is actually an hdu-list, pls rename back to hdul

Waelthus · 2023-04-27T08:45:16Z

py/picca/pk1d/prep_pk1d.py

-            flux_total_even += flexp * ivexp
-            ivar_total_even += ivexp
-            teff_even += teff_lya_exp
+    argsort = np.flip(np.argsort(np.mean(file["IV"][mask_targetid], axis=1)))


Is this line even used? I think you're redefining the var later before first usage.

Waelthus · 2023-04-27T08:54:41Z

py/picca/pk1d/prep_pk1d.py

+    ivar_mean = np.mean(file["IV"][mask_targetid][:, :], axis=1)
+    argmin_ivar = np.argmin(ivar_mean)
+    argsort = np.arange(ivar_mean.size)
+    argsort[-1], argsort[argmin_ivar] = argsort[argmin_ivar], argsort[-1]


This only moves the minimal ivar to the end (and makes the variable name confusing as it's not actually sorting the array, pls use a different name), is this a useful way to deal with things?
I.e. why not just do an actual argsort as in line 112 and use that one? Is it because the even indices are then always smaller than the odd and might there be a better way to deal with this than keeping the order of observations except for moving the worst to the end?

Yes this was some tests probably, removed

Waelthus · 2023-04-27T08:56:31Z

py/picca/pk1d/prep_pk1d.py


+    n_exp = len(flux)
+    if n_exp < 2:


shouldn't this be done at the very beginning of the routine to avoid computations

this is the first time we use the number of exp

Waelthus · 2023-04-27T08:58:41Z

py/picca/pk1d/prep_pk1d.py

+    t_even = 0
+    t_odd = 0
+    t_exp = np.sum(teff_lya)
+    for iexp in range(2 * (n_exp // 2)):


names in the previous version where closer to the picca convention of speaking variable names, e.g. index_exp instead of iexp, num_exp for nexp,...

rename in the last commit

Waelthus · 2023-04-27T09:04:59Z

py/picca/pk1d/prep_pk1d.py

+    fltotodd[w_odd] /= ivtotodd[w_odd]
+    w_even = ivtoteven > 0
+    fltoteven[w_even] /= ivtoteven[w_even]
+


please use the previous version for most things before this point as it is unchanged (except for the ivtot summation loop) and has the better variable naming convention

w_odd and w_even needed later

yes, the comment was more about notation...

Waelthus · 2023-04-27T09:12:08Z

py/picca/pk1d/prep_pk1d.py

+            fltoteven += flexp * ivexp
+            ivtoteven += ivexp
+            t_even += teff_lya_exp
+    for iexp in range(n_exp):


this can be replaced by iv_total=ivar[:n_exp].sum(axis=0). I guess similarly the loop above could be replaced by things like this:

even_inds=slice(0,2 * (n_exp // 2),2) odd_inds=slice(1,2 * (n_exp // 2),2) fltotodd=(flux[odd_inds]*ivar[odd_inds]).sum(axis=0) fltoteven=(flux[odd_inds]*ivar[even_inds]).sum(axis=0)

etc. Which would also remove all the array creation above.

Waelthus · 2023-04-27T09:13:56Z

py/picca/pk1d/prep_pk1d.py

+    if use_only_even:
+        if n_exp % 2 == 1:
+            print("Odd number of exposures discarded")
+            return None


do you want to discard odd numbers of exposures completely here? Or just discard the worst exposure? Atm you're doing the former...

we can think about those details at a later time...

Waelthus · 2023-04-27T09:16:48Z

py/picca/pk1d/prep_pk1d.py

+    elif method_alpha == "desi_time":
+        alpha = 2 * np.sqrt((t_odd * t_even) / (t_exp * (t_odd + t_even)))
+
+    diff = 0.5 * (fltoteven - fltotodd) * alpha


I'm fine with the section regarding the different methods, which of these matches the implementation in the previous version?

None, the previous version was just a test

and that test is not worth keeping even for bw-compatibility reasons, you say?

Waelthus · 2023-04-27T09:18:31Z

py/picca/pk1d/prep_pk1d.py

-                        with_correction=False,
-                        fiberid=None,
-                        log_lambda=None):
+def spectral_resolution(wdisp, with_correction=False, fiberid=None, log_lambda=None):


It would generally be nice to split functional changes from aestetic changes, I think both versions are ok here, but the old one is a little more legible...

Waelthus · 2023-06-01T08:52:07Z

@corentinravoux: Any news on this?

corentinravoux · 2023-06-01T09:50:19Z

I forgot this branch, I have make appropriate changes

Waelthus

This is ok from my side, @corentinravoux: please look at the comments regarding use_only_even and the definition of the slices.
I think it's fine to do it as is, but personally would never set the use_only_even as that discards spectra, not some last exposure, which is done automatically without the flag. If we'd actually want to use all exposures for the odd case we'd need to add functionality in a later PR.

Merging this now!

Waelthus · 2023-06-01T15:37:24Z

py/picca/pk1d/prep_pk1d.py

+    if use_only_even:
+        if n_exp % 2 == 1:
+            print("Odd number of exposures discarded")
+            return None


we can think about those details at a later time...

Waelthus · 2023-06-01T15:45:21Z

py/picca/pk1d/prep_pk1d.py

+    time_exp = np.sum(teff_lya)
+
+    even_inds = slice(0, 2 * (num_exp // 2), 2)
+    odd_inds = slice(1, 2 * (num_exp // 2), 2)


note that as written this will automatically discard the last exposure if an odd total number of exposures is there, so potentially we wouldn't need the use_only_even above at all (except that that currently does remove objects not exposures)...
If we actually wanted to include the case of odd-exposures fully, we'd need to use even_inds = slice(0, 2*(num_exp//2) + 1, 2)

Waelthus · 2023-06-01T15:47:04Z

py/picca/pk1d/prep_pk1d.py

+    fltotodd[w_odd] /= ivtotodd[w_odd]
+    w_even = ivtoteven > 0
+    fltoteven[w_even] /= ivtoteven[w_even]
+


yes, the comment was more about notation...

corentinravoux added 2 commits April 26, 2023 12:20

black formatting

b0698d2

latest diff estimation test implemented in Pk1d_DESI_changes branch

9b98518

iprafols reviewed Apr 26, 2023

View reviewed changes

Waelthus reviewed Apr 27, 2023

View reviewed changes

Waelthus assigned corentinravoux Jun 1, 2023

Waelthus and others added 2 commits June 1, 2023 10:57

Merge branch 'master' into pk1d_diff_noise_estimators

f1b5773

Michael's comments

a0d0903

formatting for pylint

a01afaa

Waelthus approved these changes Jun 1, 2023

View reviewed changes

Waelthus merged commit 0fb0717 into master Jun 1, 2023

Waelthus deleted the pk1d_diff_noise_estimators branch June 2, 2023 12:47

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Pk1d diff noise estimators #997

Pk1d diff noise estimators #997

corentinravoux commented Apr 26, 2023

iprafols left a comment

Waelthus left a comment

Waelthus Apr 27, 2023

Waelthus Apr 27, 2023

Waelthus Apr 27, 2023

corentinravoux Jun 1, 2023

Waelthus Apr 27, 2023

corentinravoux Jun 1, 2023

Waelthus Apr 27, 2023

corentinravoux Jun 1, 2023

Waelthus Apr 27, 2023

corentinravoux Jun 1, 2023

Waelthus Jun 1, 2023

Waelthus Apr 27, 2023

corentinravoux Jun 1, 2023

Waelthus Apr 27, 2023

Waelthus Jun 1, 2023

Waelthus Apr 27, 2023

corentinravoux Jun 1, 2023

Waelthus Jun 1, 2023

Waelthus Apr 27, 2023

corentinravoux Jun 1, 2023

Waelthus commented Jun 1, 2023

corentinravoux commented Jun 1, 2023

Waelthus left a comment

Waelthus Jun 1, 2023

Waelthus Jun 1, 2023

Waelthus Jun 1, 2023

Pk1d diff noise estimators #997

Pk1d diff noise estimators #997

Conversation

corentinravoux commented Apr 26, 2023

iprafols left a comment

Choose a reason for hiding this comment

Waelthus left a 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

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

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

Waelthus commented Jun 1, 2023

corentinravoux commented Jun 1, 2023

Waelthus left a 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