Adding the option to impose positive values on fitted eclipse maps #559
Conversation
|
I've fixed the conflict in S5_lightcurve_fitting/differentiable_models/StarryModel.py where the main branch was updated to allow negative rp; the positive maps branch now allows that to happen when pixel sampling isn't used |
|
Thanks for submitting this Mark - I'm excited to get this incorporated! Is it alright with you if I make a few edits directly to your PR? |
Doing this just to reduce the number of possible ECF settings and since the general intention behind the force_positivity parameter is the same between sinusoids and starry.
|
Yes, go for it! |
Adding meta.exoplanet_first (just like meta.lsq_first for emcee) to generally allow running an optimizer before sampling. Tweaking code to always save the PyMC3 trace (this will allow for resuming the sampling later if needed). Setting plotting verbosity to 1 for plot_eclipse_map figure and removing the meta.record_map setting.
Also finishing removal of meta.record_map (to always record map if pixel_ydeg is in paramtitles)
src/eureka/S5_lightcurve_fitting/differentiable_models/StarryModel.py
Outdated
Show resolved
Hide resolved
| p = pm.Normal("p", mu=0.2*amp, sd=0.2*amp, | ||
| shape=(self.npix, )) |
There was a problem hiding this comment.
Is a prior of 0.2fp ± 0.2fp really broad enough to be minimally informative and be hardcoded?
There was a problem hiding this comment.
I've changed this to be set by independent "pixel_prior_mean" and "pixel_prior_width" variables in the run parameters file; I didn't use the default free parameter setup as it doesn't include LogNormal priors so this seemed simpler. Using the peak value of the light curve (same as the fp parameter) works fine for both the mean and width parameters. I wouldn't expect a need to use anything different for either parameter, but this gives the option if needed
| if self.force_positivity: | ||
| p = pm.LogNormal("p", mu=np.log(0.3*amp), tau=1.0, | ||
| shape=(self.npix,)) |
There was a problem hiding this comment.
Similarly, is this prior really broad enough to be minimally informative and be hard coded?
There was a problem hiding this comment.
See above comment, this is now settable in the run parameters file
| # Set up amplitude differently if using pixel sampling | ||
| if 'pixel_ydeg' in self.paramtitles: | ||
| amp = temp.fp |
There was a problem hiding this comment.
Have you confirmed that the retrieved fp values are the same between sinusoids and spherical harmonics? It's important that the final retrieved value be the actual eclipse depth and not some messy normalization factor that's hard to interpret
There was a problem hiding this comment.
I've updated this to use exactly the same line in both cases, which is actually the same as before as the planetary flux from the planet_map2 defined at this stage is 1. I actually found better results fixing the fp value in the input parameters (to a sensible value), otherwise it's degenerate with the mean values of the pixels themselves. I don't know if it's possible to include the actual eclipse depth as a meaningful parameter when defining the emission in terms of pixels. Would it be a problem to do things this way (fixing fp just to give a sensible amplitude value to sample pixel values around) instead of retrieving the actual fp?
There was a problem hiding this comment.
In fact, the way it's set up at the moment, fp is irrelevant when fitting a map with pixel sampling. We could define the pixel values as normalised by an fp value which we fit for, but they're still entirely degenerate with it which slows down the sampling for no benefit. I think the best way to do it is to use the "pixel_prior_mean" as the equivalent of fp, where the user puts in the expected planetary flux magnitude, then to post-process the actual eclipse depth if it's needed
|
I found that the pixel trace saving in gradient_fitters.py wasn't working as it needed converting to arviz format before saving (did it work for you previously? It could be different dependency versions) so I changed it to convert to arviz and save. I also made it only save if it's calculating a pixel map (as the rest of the trace is already saved elsewhere in the nuts_samples.h5 file without the pixel values) |
|
I've also added the option to output the various parts of the calculated lightcurve from plot_s5.py with a "record_plot_data" switch, which I've needed for plotting things like residual eclipse mapping signals, as I don't think they're easily accessible from any of the default output options |
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #559 +/- ##
==========================================
- Coverage 58.06% 57.17% -0.89%
==========================================
Files 97 98 +1
Lines 11855 12097 +242
==========================================
+ Hits 6884 6917 +33
- Misses 4971 5180 +209 ☔ View full report in Codecov by Sentry. |
The default Starry-based eclipse mapping method fits the spherical harmonic coefficients directly, which allows maps containing regions of negative emission. In some circumstances it is desirable to require positive maps only.
This PR therefore adds the "pixel sampling" method described in https://starry.readthedocs.io/en/latest/notebooks/PixelSampling/, which is activated by adding a "pixel_ydeg" independent EPF parameter. The fitted map can then be required to be positive with a True "force_positive_map" ECF flag, which assigns a LogNormal prior to each pixel. The new "record_map" ECF flag allows the posterior distribution of the fitted maps to be recorded and plotted.
If "pixel_ydeg" is present and "force_positive_map" is False, each pixel is assigned a Normal prior allowing negative values, which gives similar results to the original method.
If this looks OK, I'm happy to add the new parameters to the documentation. I have tested fits with a variety of pixel distributions and flag combinations on my own machine.