[bump minor] Cross exposure

[corentinravoux]

Large modifications of the pk1d routines, and some modifications at the delta extraction. The main point of the PR is to implement an estimator of the one-dimensional power spectrum based on cross-correlations between different exposures of the same quasar. This PR can be extended to the computation of P3D (P_cross in particular).

For P3D: Just use picca_Pk1d.py, and the new individual pk1d file contains the deltas' Fourier transforms. For P1D cross-exposures: first use picca_delta_extraction.py with the options (use non-coadded spectra = True and keep single exposures = True) in [data], and calculate the Fourier transforms of the deltas with picca_Pk1d.py. Then, I've added a picca_Pk1D_cross_exposure.py script that estimates individual P1D cross-exposures.

The hope of this PR is to have a proper P1D cross-exposure estimator, and to start the process of implementing P3D estimators directly inside picca.

[Waelthus]

merged master back into here to enable the changes from v9, please double check... removed the "unblind" option. Not sure if that was from the starting point of this PR or one of your changes though...

[corentinravoux]

Yes this seems ok for me, I did not add or removed the unblind strategy. And the only conflict I saw was only on the addition of this option.

[corentinravoux]

@Waelthus @iprafols I would like to have your opinion especially on the modifications of delta extraction on which I worked on. Since it does not only impact the P1D.

[iprafols]

I did not check the pk1d changes in a lot of detail, let me know when you are ready for me to do so.

Also, on the changes of delta_extraction, things are starting to get more and more complicated. Maybe we should consider splitting the class DesiData into DesiData3d and DesiData1d, and have the DesiData1d inherit from the DesiData3d. I'm not saying we should do this for this PR but I do think we should discuss this at the next telecon

[corentinravoux]

Yes we could separate the reading between 3D and 1D. However, I am wondering if the changes done on the 3D and 1D analysis at this stage might not be beneficial for the other analysis? Maybe the cross-exposure could be used for the 3D one day ?

[corentinravoux]

Sorry for the long time of this PR, but I did not forgot your comments @Waelthus ! I will first try to solve the DLA masking issue, and I will come back to address those

[corentinravoux]

@Waelthus I checked your comment: I agree on the complex part, and I will implement the small changes you suggest. I am currently working (when I have time) on the proper separation between TARGETID and EXPOSUREID, it will be usefull for DLA masking and make the code cleaner.

[Waelthus]

sounds like a plan to me

[corentinravoux]

This is a working version: the spectra are read independently for different exposure, the masks are applied to individual deltas, and the continuum fitting is performed on the coadd of the exposures for each quasar.

[corentinravoux]

Coverals is giving a failing check, @Waelthus @iprafols do you think this is because I should add a test for cross-exposure ?

[Waelthus]

Limits for coveralls are set to fail when coverage decreases by > 1%. So this is the expected behaviour given a 1.2% change in coverage. The coveralls output shows that the most strongly affected files are py/picca/pk1d/compute_pk1d.py (most lines modified, but coverage actually rose) and py/picca/delta_extraction/data_catalogues/desi_data.py (mostly the verify_exposures_shape function is completely untested, but also some of the branchings due to the new options). So if you could add a test for cross-exposure configs in delta extraction, coveralls should be fine again. An actual end-to-end test of the cross-exposure pk1d calculation and post-processing is not needed for coveralls to pass, but would be useful (at the very latest when the cross exposure calculation actually gives a trustable result).

[corentinravoux]

Ok thanks! My plan to finish this PR which starts to be long:

• Implement cross-exposure test (minimal or end-to-end)
• Store the fft deltas in complex
• Test on iron

[Waelthus]

in case the current analysis routines are expected to already be able to do things properly (sry, didn't check), it might be worth changing the order (also so that the more interesting parts come first), i.e.

• initial testing the current on iron (or any smaller data set)
• implementing a test
• changing the routines to be in complex, where now you have the test to compare if results would change
• more final analyses

[corentinravoux]

I implemented the coaddition of exposures before the delta extraction. After a test on Y1 data, and even with applying multprocessing, the coaddition of a large number of data take to much time (Wall time on one perlmutter node). I added an option to not perform this step ("coadd before delta extraction"). When it is True, the delta extraction is done on the coadd and applied to individual exposures (the proper way, but the one which takes too much time on Y1) When it is False, the delta extraction is done on one of the exposure and applied to all exposure. I am testing this scheme now on Y1.

The conclusion associated to this comment is that: the coadd by rebinning routine does not seem to be optimized, it can be something that we should look into

[corentinravoux]

I changed my mind: I think we should try to avoid coadding with picca as much as possible. I will modify the cross-exposure so that the delta extraction is done independently for every exposures.

[Waelthus]

One question regarding the coadding: shouldn't one do co-adding for continuum fitting always? Or at the very least make sure that the linear quasar continuum rescaling fit to each object is the same for different spectra with same TARGETID? Else the same spectrum might have completely different coefficients for the different exposures (which I'd guess would at the very least lead to additional continuum fitting noise). But not sure why you want to avoid co-adding in delta extraction.

Updated: Sorry maybe didn't read well enough. If you're only taking the best exposure for each spectrum at continuum fitting stage I guess it's ok, but you are doing continuum fitting on lower SNR data than you'd normally have. I'm also not perfectly sure why co-adding using picca is slow. An alternative would be to change the continuum fitting stage to run its iminuit steps on collections of all the different exposures of the same TARGETID. Maybe would first need to optimize the data model for any of this to be efficient...

[corentinravoux]

One question regarding the coadding: shouldn't one do co-adding for continuum fitting always? Or at the very least make sure that the linear quasar continuum rescaling fit to each object is the same for different spectra with same TARGETID? Else the same spectrum might have completely different coefficients for the different exposures (which I'd guess would at the very least lead to additional continuum fitting noise). But not sure why you want to avoid co-adding in delta extraction

I will summarize the way we can do this in cross exposure more, as well as what I found for now:

• Continuum fitting on only one exposure for a quasar (wrong): This was the first implementation, and the simplest one. For a given target id, all the exposures had a continuum fitting, but only the last exposure considered was saved by picca.
• Continuum fitting on picca coadd: To perform this (implemented at commit 9e89216aa467838cebce6642a41d59839b66434b), picca needs to rebin all spectra individually, and separately coadd them by rebinning. It works well on a small sample. The problem is for iron, I am hitting Walltime, even with a multiprocessing of the rebinning step.
• Continuum fitting on coadd: We can read both pipeline coadd and exposures by picca but it requires to entirely change the logic of picca. And the reading time is already quite long.
• Continuum fitting on individual exposure: This is the option I was chosen. As you mentioned this might not be optimal but this is achievable

[corentinravoux]

I answered your comment before you edited it. No for now I am doing it for individual exposures separately, not the best one. This is a good idea. I can test both maybe.

[Waelthus]

The issue is that it might take a lot of changes in how I/O or cont fitting is done. Maybe it would be useful to have a look at the different continuum fitting coefficients in the run you already have for the same objects to see how much of an error is done when doing it the current way...

[corentinravoux]

After discussion, I implemented the three possible solutions for continuum fitting (all exposure independently, coadd, best exposure). I am currently testing the best exposure method on iron data.

[corentinravoux]

@iprafols @Waelthus For the coverage failure, and thus the tests of the delta extraction, I do not understand how they work. Do I need to add a test_desi_healpix_read_file for cross-exposure inside the delta_extraction/tests/data_tests.py ? Or a test_cross_exposure inside the delta_extraction/tests/expected_flux_tests.py file ? Those tests seems very unitary, is there an end-to-end delta test located in another part of picca that I did not see ?

[corentinravoux]

Additionally, with cross-exposure I need to read "spectra-" files. They are available for the tile mode but not for the healpix mode. It can be an issue for the tests I want to implement ?

[Waelthus]

Ok, I think what needs to be done to make the coverage go up is having tests analogous to add tests like the tile tests in data_tests with config settings that actually use the keep_single_exposures and PK 1D modes in a new function test_desi_tile_no_coadd or similar. I guess we'd need that also for the equivalent routines for healpix (and having the equivalent spectra files available for the healpix would be needed there, but note that the tile test should already go through most of the additional code in desi_data.py, so maybe not 100% needed to do the healpix test atm). Those should test if the data is read in correctly (including functions like verify_exposures_shape) and if you have the right number of spectra available afterwards etc...

Some functions for coadding, e.g. this are not used anywhere in the code yet, so either should be made usable via some option or be removed, I guess at survey.py as an alternative to select_best_exposure? That option needs to then have tests as well.

Things like select_best_exposures would need to have the equivalent test similar to here. In that particular case I'd probably also set the branching differently so the decision if the function needs calling is made within survey.py, while the function itself just selects the best exp. I also think that the check should not be based upon checking if a unique target list is the same as the whole target list, but instead one should carry around a variable/config option to see if those things are needed... That would clarify why things happen in a certain way without cluttering for those uniqueness tests. This should then also output results for continuum fitting iterations and thus probably as end-to-end as needed. Except for maybe the final delta computations/output-file writing

[Waelthus]

I added a skeleton for the data readin from tiles. Lets see how the coverage reacts

[Waelthus]

Ok, testing multiple ways of the tile based non-coadded analysis (i.e. with/without keeping individual spectra and one vs multiple CPUs) actually makes the coverage test pass, nevertheless would be good to have a test for the mean expected flux and healpix mode.

[corentinravoux]

I am currently testing on iron: I made one run which worked with best exposure CF, and I am running separated exposure CF to compare. Once we decide the best way to perform CF, we can implement the test on delta extraction, and maybe finally close this PR.

[corentinravoux]

Copying slack message without plots, as it is private data:

I have computed P1D on Y1 with two ways to measure the CF: CF on the best exposure, applied to all exposures CF on all exposures independently The result for independent exposures is almost perfect, no deviation to Y1 in the statistical errors bars of Y1 which are really low. And this without any modeling of the noise. With this result, I am starting wondering if it is necessary to look for correlated noise sources. The best exposure case result is not that good. There seems to be an outlier in the z=2.4 bin. But even without this one, there is a 3/4% shift for almost all scales. To understand this shift I checked the continuum of a specific quasar which have 2 independent exposures. The continuum is shown on the third plot for each exposure and for the coadd (realized in the previous runs of Y1). It seems like computing the delta of an exposure with the continuum computed on another one (either the best SNR exposure but I think also the coadd) will bias the deltas and consequently the P1D. Here, the stack of delta is not forced to zero, this solution could also help to reduce this bias, but my main conclusion is that treating exposures independently is enough.

[corentinravoux]

Since the independent exposure continuum fitting case works very well, I think I will let this solution for the PR, and remove the others solutions. @Waelthus What do you think about that ?

[Waelthus]

I think given the biases the "best exp" case might not be a good one to keep and having the 2 options of either using individual exposures for cont fitting (thus neglecting that the true QSO continuum should be ~ the same (although when there's long timelags between the obs it can also be the quasar changed a bit)) and co-adding (which I believe you've sped up). Having 2 alternative options here should be enough to allow reproducing that test later on, e.g. when writing the paper or dealing with ref comments... Agreed that the main analysis should probably be run on the indiv exp case if that's the only relatively unbiased estimate... Did you add options for selecting different exp only if they are done on different tiles as well etc? Such that teff is basically nominal for all exps (instead of having e.g. split exps as different exposures, or aren't they?)

[corentinravoux]

Ok, so for the last step of this PR, I will try to had an option which defines how the CF is performed, and it will also allows to separete the standard case to the cross-exposure case. For your last question, this is an interesting point. For now I am working on the healpix, and after checking on the spectra directory, a very large portion of the exposures are on the same tile (which is weird by the way). So, I do not think selecting on this criteria would be optimal. Another option might be to cut on the SNR scores in the spectra ? to remove exposures with EFFTIME lower than the nominal time ?

[corentinravoux]

After consultation with the desi data team, the best solution to get rid of both splitted exposures and correlated noise might be to keep only one repetition of TARGETID+NIGHT. I do not know if we add it to this PR.

[Waelthus]

What does this mean in here? Taking only one exp or rather taking a coadd of all same-TARGETID exposures in the same night, but not coadding across nights?

[corentinravoux]

I meant taking only one exposure. I think adding a coaddition step during the night might be too complicated, and I do not think we will gain much from it.

[Waelthus]

I guess it could be ok to only take the best of the split exposures in that case (I guess typical splits have most of the SNR in the first spectrum [conditions worsen while observing leading to an early stopping] OR have about equal T_eff when conditions were not great, but good enough for splits), but that would have less SNR than usual... Depending on how often exposures are actually split one could also think about removing splits entirely (that would keep the sample uniform, but we'd lose more spectra)...

[corentinravoux]

With this new commit, I have put all the options we wanted. We can focus on the testing and merge now.

[corentinravoux]

Hi @Waelthus, Thanks a lot for all those fix ! I am happy with the state of the PR, do you still have in mind unresolved issues or potential options to implement ?

[iprafols]

I am fine with merging