New PyMC3 Fitter, New Starry Astrophysical Model, and Differentiable Systematics Models

[taylorbell57]

This PR:

Resolves #265 with the addition of pymc3 as an optional dependency Works on #314 with the addition of starry as an optional dependency

In addition, this PR:

• Packages pymc3, starry, and related packages together under the pymc3 optional dependency collection
• Improves the formatting of the ecf.rst documentation page
• Adds clarity on how to install collections of optional dependencies
• Copies the input input_meta and s#_meta variables to ensure they are not changed by the next stage
• Fixes some issues with WFC3's drift_2D function and separate_direct function
• Fixes a bug with the calculation of meta.wave_low and meta.wave_hi
• Makes sure filename_S3_SpecData and filename_S4_LCData aren't written over during mergeevent
• Duplicates all non-differentiable systematic models (aside from the GP model) as differentiable models that can be used with pymc3. In the future I'd like to merge these implementations, but this was non-trivial and I didn't want to delay the addition of this important improvement
• Adds starry as the only differentiable astrophysical model. This one model handles transits, eclipses, eclipse mapping, and phase curves.
• Adds the PyMC3 No U-Turns Sampler (NUTS - called nuts in the ECF) and the differentiable gradient descent optimization as implemented in the exoplanet package (called exoplanet in the ECF). These are added in gradient_fitters.py while the old fitters remain in fitters.py. This also comes with a new Figure 5305 which shows the NUTS trace as 1D histograms and temporal walker evolution.
• Adds a new force_positivity ECF setting to Stage 5 to allow users to force positive sinusoid_pc phase curves.
• Adds a test for the exoplanet and nuts fitters using a partial phase curve fit to the synthetic MIRI eclipse observations.

Future improvements will need to include:

• Ideally merging differentiable and non-differentiable code to avoid duplicated code
• Better handling shared fits to not exceed the maximum number of fitted parameters in PyMC3. PyMC3 accepts parameters that are arrays so this shouldn't be too hard, but it'll take some work to rework the formatting from things like rp_0 to rp[0]
• Add an example figure for Figure 5305 to the documentation page

In my experience, PyMC3 NUTS fits have been as fast or significantly faster than dynesty fits, especially when not fitting for eclipse mapping signals or phase variations which are more computationally intensive. I have encountered the occasional issue where quite large datasets can't have two chains fitted in parallel and instead need to have the chains run in series, but I don't believe this could be an issue with our code.

I have successfully tested the code with and without pymc3 and related packages installed on my computer. I see no flake8 issues on my computer. All new functions have docstrings and new ECF and EPF parameters have been documented.

[codecov-commenter]

Codecov Report

Merging #492 (71cdd6b) into main (c0357d8) will increase coverage by 0.82%. The diff coverage is 64.96%.

@@            Coverage Diff             @@
##             main     #492      +/-   ##
==========================================
+ Coverage   58.81%   59.63%   +0.82%     
==========================================
  Files          78       87       +9     
  Lines        9142    10079     +937     
==========================================
+ Hits         5377     6011     +634     
- Misses       3765     4068     +303     

Impacted Files	Coverage Δ
src/eureka/S3_data_reduction/nircam.py	72.39% <0.00%> (ø)
src/eureka/S4_generate_lightcurves/generate_LD.py	12.00% <0.00%> (-0.25%)	:arrow_down:
...reka/S5_lightcurve_fitting/models/CentroidModel.py	28.57% <ø> (+0.79%)	:arrow_up:
...ka/S5_lightcurve_fitting/models/PolynomialModel.py	97.50% <ø> (-0.23%)	:arrow_down:
src/eureka/lib/citations.py	100.00% <ø> (ø)
src/eureka/lib/readECF.py	72.09% <ø> (ø)
src/eureka/S6_planet_spectra/s6_spectra.py	49.93% <8.07%> (-9.71%)	:arrow_down:
...htcurve_fitting/differentiable_models/StepModel.py	26.19% <26.19%> (ø)
src/eureka/S3_data_reduction/optspex.py	47.18% <33.33%> (-0.14%)	:arrow_down:
...rve_fitting/differentiable_models/CentroidModel.py	34.09% <34.09%> (ø)
... and 23 more

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

[taylorbell57]

Alright, your comments have been addressed, but in trying to use the code on what data we currently have for WASP-43b I found that compute_fn in Stage 6 didn't work for starry fits and then noticed that my simple tweak doesn't actually calculate the nightside hemisphere's flux but instead the anti-stellar point's flux, so I'll need to edit that function a bit more before this PR can be merged.

[taylorbell57]

Alright, as long as the tests pass here (already checked on my laptop), I think I've properly taken care of the compute_fn function for starry fits (and the compute_fp was already right because of how I setup the starry model).

It's not really clear what pc_amp and pc_offset should be for a generic spherical harmonic map, but I think I'll stick to just Y10 and Y11 for the pc_amp and pc_offset term and Y20 and Y22 for the pc_amp2 and pc_offset2 terms. I'll need to tweak the compute_amp terms as they're again calculated for the map and not the phase curve. Since I'm only using those Y10, Y11, Y20, and Y22 terms though, I'm just going to use Cowan+2008's simple but exact equation 7. I'll try to get this done tomorrow

[kevin218]

I have a more robust solution in mind for the amplitude and offset that doesn't depend on the functional form. It involves generating a series of phase curves from the parameter distributions and solving for the max, etc.

[taylorbell57]

Yeah, so that's what I did to find the nightside flux. I guess there's two different things that we're talking about though. First, I and others might want to know the amplitude and offset of the first and second order terms separately (for easier comparison with other datasets or analyses where higher order terms weren't fitted). Second, one might want to know the peak-to-trough amplitude as well as the maximum and minimum offsets of the overall combined phase curve with all the sinusoid or Ylm terms together. I guess we should probably support both - we just need clear names for all the possible options

[taylorbell57]

The other headache is that this currently only takes into account the terms that were fitted and ignores any terms set to white_fixed or fixed...

[taylorbell57]

That last commit I just made makes the behaviour of pc_amp from starry models consistent with that for batman models. I'd suggest leaving changes to both models' pc_amp calculations to a future PR. I need to locally test for a potential bug with higher order models I encountered while trying to fit WASP-43b, and after that is addressed or unreproducible I'll ping you here for your updated review.

[taylorbell57]

Alright, I figured out that I myself was just doing something wrong and that the code is fine as is. This PR is ready for your review again, @kevin218