🌈 Models¶
Often, we want to compare the data of a 🌈 to a model for the flux as a function of wavelength and time. chromatic provides a few tools to simplify performing and visualizing these comparisons. This page provides a quick tour of some of those features.
from chromatic import SimulatedRainbow, read_rainbow, version
from chromatic import plt, np, u
version()
'0.4.12'
Creating a 🌈WithModel object¶
All Rainbow objects guarantee access to 5 core quantities (wavelength, time, flux, uncertainty, ok). The RainbowWithModel object adds a model quantity (with the same shape as flux) to this list and provides new functions that make use of that model.
If we read in a Rainbow from a data file, it might not have a model set yet. We can attach an array of model values and turn a Rainbow into a RainbowWithModel object using the .attach_model() method.
data = read_rainbow("example-datasets/chromatic/ero-transit-wasp-96b.rainbow.npy")
data
<🌈(2048w, 280t)>
data_with_model = data.attach_model(
model=np.ones_like(data.flux),
planet_model=np.ones_like(data.flux),
systematics_model=np.ones_like(data.flux),
)
data_with_model
<🌈WithModel(2048w, 280t)>
The .attach_model() function requires at least the overall model be supplied; this is an array meant to represent what our flux would look like if there were no noise. We can also add additional model components, like planet_model and systematics_model in the example above, to be able to track and visualize them separately. Obviously, "ones everywhere" is probably not good a good model for a real dataset, so you'd probably want to replace the np.ones_like above with something like the outputs from a model optimization or sampling routine.
If we created a simulated dataset with a SimulatedRainbow() object, the model behind that simulation is automatically stored inside the object. Let's generate a simulated dataset and use it as an example.
simulated = SimulatedRainbow().inject_transit().inject_systematics().inject_noise()
simulated
<Simulated🌈(231w, 150t)>
Calculating 🌈 model residuals¶
The SimulatedRainbow() object inherits from the RainbowWithModel object, so it has all its powers. This includes a .residuals property, that automatically calculates flux - model based on their current values.
simulated.residuals
array([[ 0.00591736, -0.00890932, 0.00599965, ..., -0.0207381 ,
0.0100358 , -0.00560908],
[-0.00060355, -0.0114714 , -0.00806437, ..., 0.00287637,
0.00092436, -0.01388138],
[ 0.01365192, 0.0116602 , -0.00378174, ..., 0.00450487,
0.0075909 , 0.00846986],
...,
[-0.00592653, 0.00966236, -0.00917968, ..., 0.00603839,
0.01059803, 0.00052156],
[-0.01704405, -0.00325683, 0.00242677, ..., 0.00605504,
0.01062701, 0.00471661],
[-0.01342348, 0.01121836, -0.00548376, ..., 0.00457358,
-0.01019664, -0.01053731]])
In the case of our simulation, we shouldn't be surprised that the model is a good fit and the residuals look like they are drawn from a zero-mean normal distribution characterized by the uncertainty
plt.figure(figsize=(8, 2))
plt.hist((simulated.residuals / simulated.uncertainty).flatten(), bins=100)
plt.xlabel("(data - model)/uncertainty");
Visualizing 🌈 model comparisons¶
A few helpers exist for visualizing comparisons between model and data. Because they are trying to do a lot, these generally take lots of optional keyword arguments, but we've tried to make the defaults a pretty as possible.
.imshow_with_models()¶
If we have a lot of data, the most compact way to visualize it is often to show a 2D map of the data or model flux as a function of time and wavelength.
simulated.imshow_with_models();
This defaults to showing the model components systematics_model and planet_model, but we can change which models get displayed with the models= keyword argument. This won't change the residuals, which are always calculated as flux - model. The example below shows a few more common keyword arguments we often want to change.
simulated.imshow_with_models(models=["model"], xaxis="wavelength", cmap="gray");
.plot_with_model()¶
If our data are binned down to a small number of wavelengths, then it might work to plot the data to model comparison as transit light curves. The left panel below plots the data with the complete model, both of which still contain the systematics. The right panel has divided out both the flux and model arrays by the systematics model, thus making the transit much clearer.
fi, ax = plt.subplots(1, 2, sharey=True, sharex=True)
binned = simulated.bin(R=5)
binned.plot_with_model(ax=ax[0])
(binned / binned.systematics_model).plot_with_model(ax=ax[1]);
We can also plot the residuals from the model along the side.
binned.plot_with_model_and_residuals();
.animate_with_models()¶
Finally, if we want to be able to look at multiple model components with light curve plots for many wavelengths, the simplest way might be to make an animation that flips through wavelength.
binned.animate_with_models();
<IPython.core.display.Image object>
That animation function calls a helper function that you might want to use on its own, if you want to make a multicomponent lightcurve plot for a single wavelength. Here, it's being used to plot the first wavelength in the 🌈.
binned.plot_one_wavelength_with_models(i_wavelength=0);
other visualizations¶
Please remember you can always check for what other visualization or actions are available for a given 🌈 with the .help() method. All the 🌈 Visualizations methods still work for 🌈 with models attached!
If you would like another kind of a visualization that isn't shown here, please submit an Issue to discuss it.