Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Coadd1d bad spec1d #1646

Merged
merged 28 commits into from
Sep 7, 2023
Merged

Coadd1d bad spec1d #1646

merged 28 commits into from
Sep 7, 2023

Conversation

debora-pe
Copy link
Collaborator

This add a new functionality in coadd1d to reject (i.e., not coadd) spec1ds that they are bad extractions. The spec1ds are rejected based on the S/N compared to the typical S/N of all the extractions that are coadded. A new parameter has been added to guide this and is sigrej_exp.
The history in the coadd1d spectra is updated to include the info on the rejected spectra.

@debora-pe debora-pe requested a review from badpandabear August 4, 2023 23:35
nocoadd_exp = np.array([], dtype=int)

# check if there are exposures that are completely masked out, i.e., gpms = False for all pixels
masked_exp = np.where([np.all(gpm == False) for gpm in gpms])[0]
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Get rid of the where and use logical_not

bpm_exp = [np.all(np.logical_not(gpm) for gpm in gpms]

msgs.warn(f'The following {masked_exp.size} exposure(s) is/are completely masked out. '
f'It/They will not be coadded.')
[msgs.warn(f'Exposure {i}') for i in masked_exp]
# update the list of exposures that will be coadded
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See the code above, no need to use the index i or the where above, i.e. these list comprehensions can be written:

_waves = [wave for (wave, bpm) in zip(waves, bpm_exp) if not bpm_exp]

Likewise for the rest

@@ -2156,6 +2159,10 @@ def combspec(waves, fluxes, ivars, gpms, sn_smooth_npix,
Boolean mask for stacked
spectrum on wave_stack wavelength grid. True=Good.
shape=(ngrid,)
nocoadd_exp : `numpy.ndarray`_
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be changed to a list (since the inputs are list) of length next, which is True if the exposure was used and False if the exposure was not used. That is more aligned with the structure of this code.

@@ -2164,15 +2171,60 @@ def combspec(waves, fluxes, ivars, gpms, sn_smooth_npix,
_fluxes = [np.float64(flux) for flux in fluxes]
_ivars = [np.float64(ivar) for ivar in ivars]

# exposures that will not be coadded
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Change to a list. You could also call it a exp_gpm or something to make it more clear.

f'({_sigrej} sigma above the median S/N) in the stack.')
[msgs.warn(f'Exposure {i}') for i in bad_exp]
if sigrej_exp is not None:
msgs.warn('The above exposure(s) will not be coadded.')
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same feedback applies regarding the list comprehensions.

Copy link
Collaborator

@jhennawi jhennawi left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for this PR @debora-pe. I've been trying to remove masking from the core routines because it makes the control flow and the amount of code to maintain onerous if the masking is occurring in multiple places. My opinion is that whenever possible, masking should take place outside of core routines, i.e. it is understood that the core routine operates on exposures that have already been vetted. That simplifies the control flow in the core routines considerably. I've already implemented these kind of simplifications in object finding and extraction core methods. For that reason, I don't like to see masking being implemented in core routines here.

The most natural place for the vetting that you are doing is in pypeit_collate_1d (for multisite) or pypeit_flux_setup (for echelle). Another possible place for it would be in the Coadd1d class, i.e. the idea is mask in the classes, but not in core.

A counter argument would be that you need the infrastructure from this code in order to perform this vetting, but I'm not convinced about that. Basically you just need to run:

rms_sn, _ = coadd.sn_weights(fluxes, ivars, gpms, const_weights=True)

To access the rms_sn, and you need to have read in the fluxes, ivars, gpms.

I've commented on the code below, but please remove it from core.coadd.combspec and instead put in the Coadd1d multisite child or in pypeit_collate_1d.

@codecov-commenter
Copy link

codecov-commenter commented Aug 15, 2023

Codecov Report

Merging #1646 (6687db3) into develop (23c1009) will decrease coverage by 0.03%.
The diff coverage is 34.74%.

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

@@             Coverage Diff             @@
##           develop    #1646      +/-   ##
===========================================
- Coverage    40.88%   40.86%   -0.03%     
===========================================
  Files          189      189              
  Lines        43491    43575      +84     
===========================================
+ Hits         17782    17805      +23     
- Misses       25709    25770      +61     
Files Changed Coverage Δ
pypeit/coadd2d.py 8.77% <0.00%> (-0.10%) ⬇️
pypeit/core/wavecal/wvutils.py 50.00% <0.00%> (ø)
pypeit/scripts/coadd_2dspec.py 6.54% <0.00%> (ø)
pypeit/coadd1d.py 16.57% <4.08%> (-4.14%) ⬇️
pypeit/history.py 89.47% <81.39%> (-10.53%) ⬇️
pypeit/par/pypeitpar.py 96.07% <100.00%> (+<0.01%) ⬆️

📢 Have feedback on the report? Share it here.

@debora-pe
Copy link
Collaborator Author

Tests pass.

Test Summary
--------------------------------------------------------
--- PYTEST PYPEIT UNIT TESTS PASSED  234 passed, 72 warnings in 315.48s (0:05:15) ---
--- PYTEST UNIT TESTS PASSED  123 passed, 195 warnings in 1550.65s (0:25:50) ---
--- PYTEST VET TESTS PASSED  53 passed, 56 warnings in 837.68s (0:13:57) ---
--- PYPEIT DEVELOPMENT SUITE PASSED 215/215 TESTS  ---
Testing Started at 2023-08-17T14:41:38.922491
Testing Completed at 2023-08-18T00:01:05.178846
Total Time: 9:19:26.256355

pypeit/coadd1d.py Outdated Show resolved Hide resolved
# set sn_smooth_npix if not provided
if self.par['sn_smooth_npix'] is None:
# same as in coadd.multi_combspec()
sn_smooth_npix = int(np.round(0.1 * (np.sum([np.sum(wave > 1.0) for wave in _waves]) / self.nexp)))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See my feedback on your previous version of this PR. You don't need to assign sn_smooth_pix since that is used for the weights, not for the rms_sn, for which the whole spectrum is used. So all of this sn_smooth_pix should be removed. You should instead write:

rms_sn, _ = coadd.sn_weights(fluxes, ivars, gpms, const_weights=True)

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry I missed that. I removed sn_smooth_pix.
Thanks!

msgs.error('At least 2 exposures are required for coadding.')

# check if there is any bad exposure by comparing the rms_sn with the median rms_sn among all exposures
if len(_fluxes) > 2:
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This masking based on a statistical analysis of the rms_sn needs to NOT be the default behavior. My suggestion is that you set sigrej_exp to default to None. And that you only carry out the code below if sigrej_exp is not None.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Masking is not done if sigrej_exp=None. You can see below:

                if self.par['sigrej_exp'] is not None:
                    msgs.warn('The above exposure(s) will not be coadded.')

There is still some stats on the S/N, which is printed in the terminal, when sigrej_exp=None. I think it's useful to have that information, but the masking is not applied.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay I agree printing the information is useful. I would some text to the warning if sigrej_exp is not None that indicates the exposure is being masked. Right now the user won't know it is masked from watching the terminal outputs.

defaults['sn_smooth_npix'] = None
dtypes['sn_smooth_npix'] = [int, float]
descr['sn_smooth_npix'] = 'Number of pixels to median filter by when computing S/N used to decide how to scale ' \
'and weight spectra. If set to None (default), the code will determine the effective ' \
'number of good pixels per spectrum in the stack that is being co-added and use 10% of ' \
'this neff.'

defaults['sigrej_exp'] = None
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For the default of None, you should not check rms_sn and should not mask in any way.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See my response above.

Copy link
Collaborator

@jhennawi jhennawi left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @debora-pe,

Thanks. The PR look great. I have only minor comments that should be easy to address.

@debora-pe
Copy link
Collaborator Author

I re-run the tests, since I haven't done in awhile, and all is good.

Test Summary
--------------------------------------------------------
--- PYTEST PYPEIT UNIT TESTS PASSED  234 passed, 72 warnings in 313.48s (0:05:13) ---
--- PYTEST UNIT TESTS PASSED  123 passed, 195 warnings in 1510.26s (0:25:10) ---
--- PYTEST VET TESTS PASSED  53 passed, 55 warnings in 836.96s (0:13:56) ---
--- PYPEIT DEVELOPMENT SUITE PASSED 215/215 TESTS  ---
Testing Started at 2023-08-31T19:53:35.476086
Testing Completed at 2023-09-01T05:08:37.182021
Total Time: 9:15:01.705935

@debora-pe
Copy link
Collaborator Author

merging this.

@debora-pe debora-pe merged commit 0187781 into develop Sep 7, 2023
23 checks passed
@debora-pe debora-pe deleted the coadd1d_badspec1d branch September 7, 2023 23:29
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants