🌈 Actions¶
There are number of methods or actions that any Rainbow (=🌈) object can do. By learning the vocabulary of just a few of these methods, you can build up some fairly complicated stories for working with data. In general, most of these actions return a 🌈 object that has been modified in one way or another, so you can keep adding actions after actions after actions. To show how these work, we'll create some simulated 🌈 objects and try a few:
from chromatic import SimulatedRainbow, version
from chromatic import plt, np, u
plt.matplotlib.rcParams["figure.figsize"] = (8, 3)
version()
'0.4.12'
By itself, the SimulatedRainbow() object just creates an empty 🌈 with a particular wavelength and time grid. However, we can use a suite of actions that inject signals and/or noise sources to build of semi-realistic simulated datasets.
empty = SimulatedRainbow(wlim=[0.1, 5] * u.micron, tlim=[-0.1, 0.1] * u.day)
🌈.inject_transit()¶
It's useful to have a way to inject a simulated transit of an imaginary exoplanet. If you're trying out some new fitting algorithm, you can inject a transit with known properties and make sure you recover it.
with_transit = empty.inject_transit()
with_transit.imshow();
You can customize some of the model parameters. See the docstring for more details, but the most basic would be to provide an array of planet radii corresponding to each wavelength.
with_transit = empty.inject_transit(planet_radius=np.linspace(0.2, 0.1, empty.nwave))
with_transit.imshow();
The parameters used for the transit model will be stored in .metadata['injected_transit_parameters'].
🌈.inject_spectrum()¶
This function injects a static stellar spectrum into all the fluxes for a 🌈, using the Husser et al. (2013) PHOENIX model grid. It will try to automatically download all the files you need, when you need them.
with_spectrum = with_transit.inject_spectrum(
temperature=5800 * u.K, logg=4.43, metallicity=0.0
)
with_spectrum.imshow();
🌈🤖 None of the pre-existing flux values were 1, which hints at the possibility that there might already be a spectrum in them. Please watch out for weird units or values!
🌈.inject_systematics()¶
This function injects imaginary systematic noise sources into a 🌈. High-precision observations of transiting exoplanets often encounter systematics that are correlated in complicated ways in time or wavelength or with various hidden quantities. This function injects a linear combination of real or imagined time-like, wave-like, and/or flux-like quantities (see its docstring for details).
with_systematics = with_spectrum.inject_systematics()
with_systematics.normalize().imshow();
The parameters and equations used for the systematics model will be stored in .metadata['systematics_components'] and .metadata['systematics_equation']. The independent variables needed for these equations are stored in the 🌈, so it should be possible to perfectly recreate the systematic model.
🌈.inject_noise()¶
We can inject noise, for example to simulate the observing the same system at a greater distance or with a smaller telescope. One way to set the noise level is with the signal_to_noise keyword argument.
with_noise = with_systematics.inject_noise(signal_to_noise=100)
with_noise.normalize().imshow();
Another way to set the noise level is to specify a typical number of photons expected per bin with the number_of_photons keyword argument.
with_noise = with_systematics.normalize().inject_noise(number_of_photons=1e4)
with_noise.imshow();
Since each 🌈 action returns a 🌈, multiple actions can be linked togther into a single command.
r = (
SimulatedRainbow()
.inject_transit()
.inject_spectrum()
.inject_systematics()
.inject_noise()
.normalize()
)
🌈🤖 None of the pre-existing flux values were 1, which hints at the possibility that there might already be a spectrum in them. Please watch out for weird units or values!
🌈.bin()¶
While we should generally try to avoid fitting to binned data when possible, there will often be times where it's helpful to bin to particular grid of wavelengths and/or times. You can do this using the .bin() function, which we'll apply here to the simulated 🌈 we just created.
def summarize(x):
print(
f"""
{x} is a {type(x).__name__}.
It has a {x.nwave} wavelengths and {x.ntime} times.
Its 5 first wavelengths:{np.round(x.wavelength[:5], 2)}
Its 5 first times:{np.round(x.time[:5].to('hour'), 2)}
"""
)
summarize(r)
<Simulated🌈(231w, 150t)> is a SimulatedRainbow.
It has a 231 wavelengths and 150 times.
Its 5 first wavelengths:[0.5 0.51 0.51 0.52 0.52] micron
Its 5 first times:[-2.5 -2.47 -2.43 -2.4 -2.37] h
To bin in wavelength, it can take the following inputs:
dw=to bin in wavelength to a particular $d\lambda$ width. This will create a linear grid in wavelength.R=to bin in wavelength to particular $R = \lambda/d\lambda$. This will create a logarithmic grid in wavelength.wavelength=to bin to any custom wavelength grid specified by its centers; the edges will be guessed from the spacing between these edges.wavelength_edges=to bin to any custom wavelength grid specified by its edges; the centers will be guessed as the midpoints between these edges. (The number of binned wavelengths will be 1 fewer than the number of edges.)nwavelengths=to bin by a fixed number of adjacent wavelengths (as in "bin every N wavelengths together"), starting from the first wavelength.
b = r.bin(dw=0.5 * u.micron)
summarize(b)
<Simulated🌈(10w, 150t)> is a SimulatedRainbow.
It has a 10 wavelengths and 150 times.
Its 5 first wavelengths:[0.5 1. 1.5 2. 2.5] micron
Its 5 first times:[-2.5 -2.47 -2.43 -2.4 -2.37] h
b = r.bin(R=10)
summarize(b)
<Simulated🌈(24w, 150t)> is a SimulatedRainbow.
It has a 24 wavelengths and 150 times.
Its 5 first wavelengths:[0.5 0.55 0.61 0.67 0.75] micron
Its 5 first times:[-2.5 -2.47 -2.43 -2.4 -2.37] h
b = r.bin(wavelength=np.linspace(1, 2, 6) * u.micron)
summarize(b)
<Simulated🌈(6w, 150t)> is a SimulatedRainbow.
It has a 6 wavelengths and 150 times.
Its 5 first wavelengths:[1. 1.2 1.4 1.6 1.8] micron
Its 5 first times:[-2.5 -2.47 -2.43 -2.4 -2.37] h
b = r.bin(wavelength_edges=np.linspace(1, 2, 6) * u.micron)
summarize(b)
<Simulated🌈(5w, 150t)> is a SimulatedRainbow.
It has a 5 wavelengths and 150 times.
Its 5 first wavelengths:[1.1 1.3 1.5 1.7 1.9] micron
Its 5 first times:[-2.5 -2.47 -2.43 -2.4 -2.37] h
b = r.bin(nwavelengths=10)
summarize(b)
<Simulated🌈(23w, 150t)> is a SimulatedRainbow.
It has a 23 wavelengths and 150 times.
Its 5 first wavelengths:[0.52 0.58 0.64 0.71 0.78] micron
Its 5 first times:[-2.5 -2.47 -2.43 -2.4 -2.37] h
To bin in time, it can take the following inputs:
dt=to bin in time to a particular $dt$ width. This will create a linear grid in time.time=to bin to any custom time grid specified by its centers; the edges will be guessed from the spacing between these edges.time_edges=to bin to any custom time grid specified by its edges; the centers will be guessed as the midpoints between these edges. (The number of binned times will be 1 fewer than the number of edges.)ntimes=to bin by a fixed number of adjacent times (as in "bin every N times together"), starting from the first time.
b = r.bin(dt=0.25 * u.hour)
summarize(b)
<Simulated🌈(231w, 21t)> is a SimulatedRainbow.
It has a 231 wavelengths and 21 times.
Its 5 first wavelengths:[0.5 0.51 0.51 0.52 0.52] micron
Its 5 first times:[-2.5 -2.25 -2. -1.75 -1.5 ] h
b = r.bin(time=np.linspace(-1, 1, 6) * u.hour)
summarize(b)
<Simulated🌈(231w, 6t)> is a SimulatedRainbow.
It has a 231 wavelengths and 6 times.
Its 5 first wavelengths:[0.5 0.51 0.51 0.52 0.52] micron
Its 5 first times:[-1. -0.6 -0.2 0.2 0.6] h
b = r.bin(time_edges=np.linspace(-1, 1, 6) * u.hour)
summarize(b)
<Simulated🌈(231w, 5t)> is a SimulatedRainbow.
It has a 231 wavelengths and 5 times.
Its 5 first wavelengths:[0.5 0.51 0.51 0.52 0.52] micron
Its 5 first times:[-0.8 -0.4 0. 0.4 0.8] h
b = r.bin(ntimes=10)
summarize(b)
<Simulated🌈(231w, 15t)> is a SimulatedRainbow.
It has a 231 wavelengths and 15 times.
Its 5 first wavelengths:[0.5 0.51 0.51 0.52 0.52] micron
Its 5 first times:[-2.35 -2.02 -1.68 -1.35 -1.02] h
You can combine to bin in both wavelength and time in the same step.
b = r.bin(dw=100 * u.nm, dt=10 * u.minute)
summarize(b)
<Simulated🌈(46w, 31t)> is a SimulatedRainbow.
It has a 46 wavelengths and 31 times.
Its 5 first wavelengths:[0.5 0.6 0.7 0.8 0.9] micron
Its 5 first times:[-2.5 -2.33 -2.17 -2. -1.83] h
fi, ax = plt.subplots(1, 2, sharex=True, figsize=(8, 3), constrained_layout=True)
r.imshow(ax=ax[0], vmin=0.98, vmax=1.02)
plt.title("Unbinned")
plt.xlabel("")
b.imshow(ax=ax[1], vmin=0.98, vmax=1.02)
plt.title("Binned");
🌈.fold()¶
If our times are in units of BJD, it might be helpful to phase-fold them to a particular transit period and epoch. This function is a small wrapper to convert time to being measured relative to the center of some periodic event (like a transit).
from astropy.time import Time
r.time += Time.now().jd * u.day
r.time
r.fold(period=1.234 * u.day, t0=Time.now().jd * u.day).time
🌈.get_average_lightcurve_as_rainbow()¶
Binning all wavelengths together produces an uncertainty-weighted average light curve.
r.get_average_lightcurve_as_rainbow()
<Simulated🌈(1w, 150t)>
🌈.get_average_spectrum_as_rainbow()¶
Binning all times together produces an uncertainty-weighted average spectrum.
r.get_average_spectrum_as_rainbow()
<Simulated🌈(231w, 1t)>
🌈.normalize()¶
If we're starting in units of photons detected at our telescope per wavelength, we may want to normalize a 🌈 by dividing through by its median spectrum. That's effectively all that the .normalize() action does.
r.normalize().imshow();
In some cases, you might also be curious to divide by the median light curve, to look for small variations away from the overall transit shape. You can customize whether you want to normalize in wavelength or time by suppyling the axis= keyword when you call the .normalize() function, which defaults to axis='wavelength'.
r.normalize(axis="time").imshow();
(In this case, there were no transit depth or limb-darkening variations injected into the simulated 🌈, so normalizing through time entirely erases any hint of the transit.)
🌈.remove_trends()¶
Often we need a quick way to remove smooth trends from a 🌈.
fi, ax = plt.subplots(2, 2, sharey=True, figsize=(8, 5), constrained_layout=True)
r.imshow(ax=ax[0, 0])
plt.title("original")
r.remove_trends(method="median_filter", size=(11, 5)).imshow(ax=ax[0, 1])
plt.title("median-filtered")
r.remove_trends(method="savgol_filter", window_length=11, polyorder=2).imshow(
ax=ax[1, 0]
)
plt.title("savgol-filtered")
r.remove_trends(method="differences").imshow(ax=ax[1, 1])
plt.title("first-differences");
🌈.shift()¶
This function Doppler shifts all the wavelengths in a 🌈 by a given velocity.
unshifted = SimulatedRainbow(wlim=[470, 560] * u.nm, R=5000).inject_spectrum()
shifted = unshifted.shift(300 * u.km / u.s)
for x, label in zip([unshifted, shifted], ["unshifted", "shifted"]):
plt.plot(x.wavelength, x.flux[:, 0], label=label)
plt.xlabel(f"Wavelength ({unshifted.wavelength.unit.to_string('latex_inline')})")
plt.ylabel(f"Flux ({unshifted.flux.unit.to_string('latex_inline')})")
plt.legend(frameon=False);
🌈.trim()¶
Often datasets may have blocks of wavelengths or times that are entirely bad. This function trims bad wavelengths or times off the edges of a 🌈.
r.wavelike["ok"] = np.random.uniform(0, 1, r.nwave) < 0.9
r.wavelike["ok"][:20] = False
fi, ax = plt.subplots(1, 2, figsize=(8, 3), sharey=True)
r.trim().imshow(ax=ax[1])
r.imshow(ax=ax[0]);
🌈 + 🌈¶
We can perform mathematical operations on 🌈 objects. This will apply the requested mathematical operation between the flux and model arrays, and try its best to figure out what to do with the uncertainty.
a = SimulatedRainbow().inject_noise()
b = SimulatedRainbow().inject_noise()
(a + b).imshow();
(a - b)
<Simulated🌈(231w, 150t)>
(a * b)
<Simulated🌈(231w, 150t)>
(a / b)
<Simulated🌈(231w, 150t)>
a * 2 + b / 4 - 0.1 * a.model
<Simulated🌈(231w, 150t)>
🌈[:,:]¶
We can trim a Rainbow in wavelength (first dimension) and/or time (second dimension) by slicing it in a similar way that you might to any other 2D array.
# the original array
summarize(r)
<Simulated🌈(231w, 150t)> is a SimulatedRainbow.
It has a 231 wavelengths and 150 times.
Its 5 first wavelengths:[0.5 0.51 0.51 0.52 0.52] micron
Its 5 first times:[59053349.96 59053350. 59053350.03 59053350.06 59053350.1 ] h
# using slices
summarize(r[5:10, ::2])
<Simulated🌈(5w, 75t)> is a SimulatedRainbow.
It has a 5 wavelengths and 75 times.
Its 5 first wavelengths:[0.53 0.53 0.54 0.54 0.55] micron
Its 5 first times:[59053349.96 59053350.03 59053350.1 59053350.16 59053350.23] h
# using indices
i_wavelengths = np.arange(5, 10)
i_times = np.arange(0, r.ntime, 2)
summarize(r[i_wavelengths, i_times])
<Simulated🌈(5w, 75t)> is a SimulatedRainbow.
It has a 5 wavelengths and 75 times.
Its 5 first wavelengths:[0.53 0.53 0.54 0.54 0.55] micron
Its 5 first times:[59053349.96 59053350.03 59053350.1 59053350.16 59053350.23] h
# using boolean masks
ok_wavelengths = r.wavelength < 2 * u.micron
ok_times = np.abs(r.time) < 1 * u.hour
summarize(r[ok_wavelengths, ok_times])
<Simulated🌈(139w, 0t)> is a SimulatedRainbow.
It has a 139 wavelengths and 0 times.
Its 5 first wavelengths:[0.5 0.51 0.51 0.52 0.52] micron
Its 5 first times:[] h
Viewing a 🌈's History¶
Most actions that return Rainbow objects be recorded in that object's metadata['history'] entry. This is meant to be an approximate summary of the steps that led up to the creation of the current 🌈. You can view this history by calling the .history() method.
x = (
SimulatedRainbow()
.inject_noise()
.inject_transit()
.bin(R=5, dt=5 * u.minute)
.normalize()
)
h = x.history()
print(h)
(
SimulatedRainbow(
tlim=u.Quantity(np.array([-2.5, 2.5]))*u.Unit('h'),
dt=u.Quantity(2.0)*u.Unit('min'),
wlim=u.Quantity(np.array([0.5, 5. ]))*u.Unit('micron'),
R=100)
.inject_noise(
signal_to_noise=100)
.inject_transit(
planet_radius=0.1,
method='exoplanet',
**{})
.bin_in_time(
dt=u.Quantity(5.0)*u.Unit('min'),
minimum_acceptable_ok=1,
trim=True)
.trim_times(
just_edges=True,
when_to_give_up=1,
minimum_acceptable_ok=1)
.bin_in_wavelength(
R=5,
minimum_acceptable_ok=1,
trim=True,
starting_wavelengths='1D')
.trim_wavelengths(
just_edges=True,
when_to_give_up=1,
minimum_acceptable_ok=1)
.normalize(
axis='wavelength',
percentile=50)
)
By default, it gives an almost copy-paste-able string version of the history of the Rainbow. If you look closely, you'll see that .bin has gotten split out into four steps (.bin_in_time, .trim_times, .bin_in_wavelength, .trim_wavelengths). In many cases, you should be able to approximately reproduce the actions that have gone into a Rainbow by just copying, pasting, and running the set of commands returned by .history() (or running eval to evaluate a string as Python commands).
eval(h)
<Simulated🌈(13w, 61t)>