-
Notifications
You must be signed in to change notification settings - Fork 52
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
POET transit and eclipse models, SPAM LD #609
Conversation
It would seem much cleaner to me to generalize the transit fitting code to allow multiple planets rather than requiring a different type of transit model for each transit. I have something in trappist1_phasecurves that does this in a very hacky way for the TRAPPIST-1 b+c phase curve (which has transits of b, c, and g; eclipses of b and c; and phase variations of b and c), but that code is currently very hard-coded for that system and not generalized at all (and only works for the starry code). That's not to say that adding the POET transit model isn't worthwhile (especially given the speed improvement for the 4-parameter LD law), just that doing so for the purpose of fitting two transits doesn't seem like a good tactic to me (as we're going to need to support >2 planets anyway). Additionally, it will likely be confusing to users to have the names of the transit model parameters vary between the batman and the POET implementations; if someone wants to switch their model, they'll have to rename all of the variables in their EPFs which would be annoying. And folks writing an EPF from scratch might use some weird mixture of both parameter sets by accident. Additionally, rather than adding a different phase curve model for the POET transit, I would suggest generalizing the current phase curve model to work with either the batman or the POET transit/eclipse models; the two phase curve models are mathematically the same model just coded in two different ways, so it'll save us a fair bit of duplicated code and some maintenance effort if we just merge the two functions. |
I couldn't find a reasonable way to implement two-planet fits with BATMAN, but when I get back from vacation I'll take a look at your branch to see what you did. I still want to implement POET transit and eclipse models. If we can get two-planet fitting up and running quickly, then I'm happy to adopt the same parameter names between POET and BATMAN. I'm on a bit of a time crunch to solve this problem because we're also working on a dataset with multiple planet transits. I'll take a second look at broadening the phase curve model to allow for POET and BATMAN models. That was the major concern I had. |
Everything looks good to me and all the tests are passing which is great. I'll just try doing a quick optimization fit to the TRAPPIST-1 b+c phase curve data to ensure it runs properly and check if its generally consistent with the batman model |
FYI, I've identified a couple of issues (some of which are related to this PR, some of which aren't). I'll try to separate and specify the issues related to this PR and may potentially just submit patches myself for several of them (likely as a PR that merges into this PR so you can look over the changes). One example of a current issue with this PR is that the alias between |
Old bugs include issues with multwhite fits, masked values, and nbin_plot dealing with nbin_plot=len(time) in the presence of masked values. New POET bugs include improperly handling multi-channel (or multwhite) fits (because only the values for channel 0 were ever used), and a broken alias between fp and fpfs.
gp.lnlikelihood was not recieving the residuals, so the george GP wasn't going to work properly. george GP was also not receiving only the unmasked kernel_inputs.
Older versions of numpy that are needed for pymc3 do not have np.ma.ones_like or np.ma.zeros_like, so reworking the bits of code that used those functions
Added testing for centroid models in PyMC3 and numpy codes, added S5 testing for photometry. Also pytest>=8 doesn't seem to work at all if the pymc3 code is installed, so just setting an upper-limit of <8 since we don't need newer versions of that package
Fixing new POET bugs, and old bugs with multwhite fits, masked values in fits, nbin_plot in S5, and S6 model loading
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Everything seems to work well now!
The one caveat is that fitting the TRAPPIST-1 phase curve while fitting for ecc
and w
instead of t_secondary
, the POET fit settles in a fairly different regime because it allows ecc
to get larger than the batman model does given the impact on the transit/eclipse shapes (as best as I can tell), and the POET fit also ends up with a larger GP component to soak up the residual signal. That doesn't seem to be indicative of a bug with the code, just the result of the assumption that eccentricity doesn't impact transit/eclipse shapes.
But, when I use the approach I intend to use in practice of fitting for t_secondary
instead while fixing ecc=0
(since its close to that anyway for the inner TRAPPIST-1 planets), then the POET and batman models agree quite well with the POET models being slower (247 minutes vs 140 minutes with batman) but more easily interpretable
This PR adds an alternate option to BATMAN when fitting planet transits and eclipses. This also allows users to fit two planet transits, one with BATMAN and the other with POET. The sinusoidal phase curve model also depends on using BATMAN, so I've ported over the POET phase curve model.
Tests show consistent results between BATMAN and POET LC fits. In most cases, they are similar in speed; however, POET sees a ~50x improvement in speed using when using the
4-parameter
LD model because it adopts the small-planet approximation (which is valid for RpRs < 0.1).Because I need it for the same project, I've added the ability at Stage 4 to load in limb-darkening files (columns are wavelength, u1, u2, etc.) generated by SPAM. The file can be an arbitrarily high resolution and the new code will compute a simple average for each spectroscopic channel.
What's left to do...
What I'm leaving for another day...
jit