kevin218 / Eureka

Eureka! is a data reduction and analysis pipeline intended for time-series observations with JWST.
https://eurekadocs.readthedocs.io/
MIT License
58 stars 45 forks source link

(1) Allows negative Rp/Rs for shallow transits, (2) nansums in optspec extraction to address masking issues, (3) allows for skipping/fitting only certain bins in S5, and (4) updates aperture to be centered properly #533

Closed erinmmay closed 1 year ago

erinmmay commented 1 year ago

Lots of tiny changes that I'm finally looking to merge into the main branch.

1. Allows negative Rp/Rs for shallow transits For shallow transits, we shouldn't reject steps that try negative Rp/Rs. This results in chopped-off posteriors. I simply removed that condition from the "physicality" check.

2. nansums in optspec extraction to address masking issues In the spectral extraction, the masking is setting masked values to nan. Optspec.py uses regular sums, so you end up with fully rejected columns if there was a single bad pixel/outlier, which can get out of hand if you have may bad/trimmed pixels. Switching to nansum fixes this.

3. allows for skipping/fitting only certain bins in S5 Simple change to allow input in S5 that lets you only fit certain bins, or skip certain bins. Also wraps the fitting in a try/except statement so a failed fit doesn't crash the full run.

4. updates aperture to be centered properly Based on discussions with Kevin, aperture extraction was from [pos-aperture: pos+aperture], which means you get one fewer pixel above the center than below. Updated to be centered, results in an odd number of pixels within the aperture centered at the center of the trace. So the total number of pixels is the 2*ap_hw + 1 now. Similar update for the background region. This also addresses a half pixel offset that is used for computing the center of mass for computing how to roll columns when straightening the trace which isn't removed prior to calculating the number of pixels to roll by.

codecov-commenter commented 1 year ago

Codecov Report

Merging #533 (1fb3fe3) into main (3c10926) will not change coverage. The diff coverage is 100.00%.

:mega: This organization is not using Codecov’s GitHub App Integration. We recommend you install it so Codecov can continue to function properly for your repositories. Learn more

@@           Coverage Diff           @@
##             main     #533   +/-   ##
=======================================
  Coverage   57.53%   57.53%           
=======================================
  Files          92       92           
  Lines       11076    11076           
=======================================
  Hits         6373     6373           
  Misses       4703     4703           
Impacted Files Coverage Δ
src/eureka/S3_data_reduction/nircam.py 74.86% <100.00%> (ø)
src/eureka/S3_data_reduction/optspex.py 49.05% <100.00%> (ø)
src/eureka/S3_data_reduction/plots_s3.py 68.47% <100.00%> (ø)
src/eureka/S3_data_reduction/s3_reduce.py 87.41% <100.00%> (ø)
src/eureka/S3_data_reduction/straighten.py 100.00% <100.00%> (ø)
...ureka/S5_lightcurve_fitting/models/BatmanModels.py 70.21% <100.00%> (ø)
src/eureka/S5_lightcurve_fitting/s5_fit.py 68.71% <100.00%> (ø)

:mega: We’re building smart automated test selection to slash your CI/CD build times. Learn more

taylorbell57 commented 1 year ago

Another thing to consider is what impact all of this will have on the Stage 6 code - my guess is that many of the changes in this PR are likely to break much of the Stage 6 functionality.

taylorbell57 commented 1 year ago

So, your point 2 related to Stage 3 optspex.py nan masking or np.nan functions, the main underlying issue is that the mask array is ignored in some key pieces of the standard aperture extraction which causes a cascading set of issues. While the np.nan functions you switched to here avoid entirely nan columns in the 2D lightcurves, they only partially resolve the underlying issue that masked pixels are not masked during this part of the processing. As a result, I think the solution from my PR #549 (while imperfect and still in need of some fine tuning before merging) is a better solution

taylorbell57 commented 1 year ago

@erinmmay, would you like help in getting this PR merged? I just noticed there haven't been any changes to the PR in the past 3 months, and I'm happy to help contribute some edits to get this ready for merging if you want

erinmmay commented 1 year ago

Yeah, sorry about the delay! I think the easiest/quickest thing will be to pull out the stage 5 changes for now since that requires some extra work and is the least important of these changes. I think then it's just looking at negative transits in starry (which I've never used and would definitely appreciate you looking into if you have used it) and stripping back the nansums so your changes to the masking can be done instead. If that plan sounds good I can make this a priority.

taylorbell57 commented 1 year ago

Yeah, punting your point 3 to another PR to be able to quickly merge points 1 and 4 sounds like a good idea to me. That's not to say that point 3 isn't worth pursuing, it's just going to be a bit finicky to get coded just right.

I'll quickly test a negative radius in starry to see what happens. Even if starry doesn't allow it (which it may since I don't know if all that fancy math would break with a negative radius), we can likely still make it work by manually inverting the transit signal

taylorbell57 commented 1 year ago

Alright, so while starry doesn't technically crash with negative radii, it doesn't end up giving the expected behaviour. I suggest the following which should have fairly minimial impacts on speed (you're welcome to implement it, or I can too). I'll note that I haven't tested this fully implemented or anything, but I have tested it in a Jupyter notebook and think it should work.


In eureka.S5_lightcurve_fitting.differentiable_models.StarryModel.eval() in place of the current line which reads:

lcpiece = systems[chan].flux(time)

we need to instead put:

fstar, fp = systems[chan].flux(time, total=False)
if lessthan(planet.r, 0):
    fstar = 2-fstar
    fp *= -1
lcpiece = fstar+fp

where we'll need to define greater above in the if eval: section as:

        if eval:
            lib = np
            lessthan = np.less
            systems = self.fit.systems
        else:
            lib = tt
            lessthan = tt.lt
            systems = self.systems



And in eureka.S5_lightcurve_fitting.differentiable_models.StarryModel.setup() where we have:

r=temp.rp*temp.Rs

we instead need

r=tt.abs_(temp.rp)*temp.Rs

and in eureka.S5_lightcurve_fitting.differentiable_models.StarryModel.update() where we have:

r=temp.rp*temp.Rs

we instead need

r=np.abs(temp.rp)*temp.Rs



Finally, in eureka.S5_lightcurve_fitting.differentiable_models.PyMC3Models.CompositePyMC3Model.__init__(), where we currently have:

                                elif parname in ['rp', 'per', 'scatter_mult',
                                                 'scatter_ppm', 'c0', 'r1',
                                                 'r4', 'r7', 'r10']:

we'll need

                                elif parname in ['per', 'scatter_mult',
                                                 'scatter_ppm', 'c0', 'r1',
                                                 'r4', 'r7', 'r10']:
erinmmay commented 1 year ago

@taylorbell57 I think I've successfully removed the necessary commits. double checking things right now

I've also added a new commit that addresses a plotting issue where the y-axis wasn't correct in the residual background plot. The aperture markings on that make much more sense now.

erinmmay commented 1 year ago

ahh, I see it deleted s5_fit instead of ignoring the file... can you reject that change on your end or should I add it back in?

taylorbell57 commented 1 year ago

No, I can't accept/reject just parts of a PR - it's all or nothing on the reviewer's end

erinmmay commented 1 year ago

👍 wasn't sure if you could do it in the "resolve conflicts" stage. I'll get that file added back in.

taylorbell57 commented 1 year ago

I'll add the PyMC3 edits to this PR now if that's alright with you?

erinmmay commented 1 year ago

Yep, sounds good!

taylorbell57 commented 1 year ago

@erinmmay I don't think you've set it so that reviewers can make edits to your PR - I've made edits to PRs from forks before and can't seem to do so for yours, so my assumption is that you didn't check that box. Can you see if you're able to check it now (or confirm that it's already checked and I'm just missing something on my end)?

taylorbell57 commented 1 year ago

Aha, never mind - I am able to make commits. It's just been a long time since I did that and couldn't remember how

taylorbell57 commented 1 year ago

Alright, I think this PR looks ready to me - do you want to weigh in @kevin218, or should I approve it?

taylorbell57 commented 1 year ago

I can also confirm that all of the tests pass locally as well (GitHub automated testing is broken for all our PRs right now)

kevin218 commented 1 year ago

@taylorbell57 You can go ahead and approve.