Multiband parametric fit: RAINBOW model¶
RainbowFit fits a physical blackbody + bolometric-light-curve model to multi-band
photometric time series. It is designed for extragalactic transients such as supernovae,
kilonovae, or tidal disruption events, where the SED is well approximated by a Planck
function whose temperature and luminosity evolve with time.
The method is described in Russeil et al. 2024.
Prerequisites¶
RainbowFit requires iminuit ≥ 2.21.
Install it together with all optional dependencies or separately.
We would also need matplotlib for plotting, but of course it is not needed for the fitting itself.
In [1]:
Copied!
# %pip install "light-curve[full]" matplotlib
# OR
# %pip install light-curve iminuit matplotlib
# %pip install "light-curve[full]" matplotlib
# OR
# %pip install light-curve iminuit matplotlib
In [2]:
Copied!
import numpy as np
import matplotlib.pyplot as plt
from light_curve import RainbowFit
import numpy as np
import matplotlib.pyplot as plt
from light_curve import RainbowFit
Generating synthetic supernova data¶
We simulate a Bazin-shaped bolometric light curve with a sigmoid temperature evolution,
observed in four LSST-like bands (g r i z). Observations are flux densities — not
magnitudes — because RainbowFit works in linear flux space.
The sigmoid temperature is parametrized by the mid temperature T = $(T_\text{min}+T_\text{max})/2$
and a dimensionless relative amplitude T_amplitude = $(T_\text{max}-T_\text{min})/(T_\text{max}+T_\text{min})$,
so $T_\text{max} = T(1 + \texttt{T\_amplitude})$ and $T_\text{min} = T(1 - \texttt{T\_amplitude})$. Here the
source cools between 15 000 K and 5 000 K, i.e. T = 10 000 K and T_amplitude = 0.5.
Physical parameters¶
| Parameter | Symbol | Value |
|---|---|---|
| Peak time | $t_0$ | 60 000 (MJD) |
| Amplitude | $A$ | 1.0 |
| Rise time | $t_\text{rise}$ | 5 d |
| Fall time | $t_\text{fall}$ | 30 d |
| Mid temperature | $T = (T_\text{min}+T_\text{max})/2$ | 10 000 K |
| Relative amplitude | $T_\text{amplitude} = \frac{T_\text{max}-T_\text{min}}{T_\text{max}+T_\text{min}}$ | 0.5 (i.e. $T_\text{min}$ = 5 000 K, $T_\text{max}$ = 15 000 K) |
| Colour timescale | $t_\text{color}$ | 10 d |
In [3]:
Copied!
rng = np.random.default_rng(0)
# Effective wavelengths in Angstrom for griz
band_wave_aa = {"g": 4770.0, "r": 6231.0, "i": 7625.0, "z": 9134.0}
# True parameters
reference_time = 60_000.0
amplitude = 1.0
rise_time = 5.0
fall_time = 30.0
Tmax = 15_000.0 # hot (early) temperature
Tmin = 5_000.0 # cool (late) floor
t_color = 10.0
# Sigmoid temperature fit parameters: mid temperature and relative amplitude
T_mid = 0.5 * (Tmin + Tmax) # = 10 000 K
T_amplitude = (Tmax - Tmin) / (Tmax + Tmin) # = 0.5
# Per-band baseline flux (constant offset)
baselines = {b: 0.3 * amplitude + rng.exponential(scale=0.3 * amplitude) for b in band_wave_aa}
# Build the feature object first so we can use .model() to generate data
feature = RainbowFit.from_angstrom(
band_wave_aa,
with_baseline=True,
bolometric="bazin",
temperature="sigmoid",
)
# True parameter vector, in feature.names order:
# common + bolometric + temperature + baselines + chi2_placeholder
true_params = [
reference_time, amplitude, rise_time, fall_time, # Bazin bolometric
T_mid, T_amplitude, t_color, # sigmoid temperature: T (mid), T_amplitude, t_color
*baselines.values(), # per-band baselines
1.0, # reduced chi^2 placeholder
]
# Random observation times covering ±3 rise/fall windows
n_obs = 1000
t = np.sort(rng.uniform(reference_time - 3 * rise_time, reference_time + 3 * fall_time, n_obs))
band = rng.choice(list(band_wave_aa), size=n_obs)
# True model flux
flux = feature.model(t, band, *true_params)
# Poisson-like noise: S/N=10 at minimum flux
flux_err = np.sqrt(flux * np.min(flux)) / 10.0
flux = flux + rng.normal(0.0, flux_err)
print(f"Generated {n_obs} observations in bands {list(band_wave_aa)}")
rng = np.random.default_rng(0)
# Effective wavelengths in Angstrom for griz
band_wave_aa = {"g": 4770.0, "r": 6231.0, "i": 7625.0, "z": 9134.0}
# True parameters
reference_time = 60_000.0
amplitude = 1.0
rise_time = 5.0
fall_time = 30.0
Tmax = 15_000.0 # hot (early) temperature
Tmin = 5_000.0 # cool (late) floor
t_color = 10.0
# Sigmoid temperature fit parameters: mid temperature and relative amplitude
T_mid = 0.5 * (Tmin + Tmax) # = 10 000 K
T_amplitude = (Tmax - Tmin) / (Tmax + Tmin) # = 0.5
# Per-band baseline flux (constant offset)
baselines = {b: 0.3 * amplitude + rng.exponential(scale=0.3 * amplitude) for b in band_wave_aa}
# Build the feature object first so we can use .model() to generate data
feature = RainbowFit.from_angstrom(
band_wave_aa,
with_baseline=True,
bolometric="bazin",
temperature="sigmoid",
)
# True parameter vector, in feature.names order:
# common + bolometric + temperature + baselines + chi2_placeholder
true_params = [
reference_time, amplitude, rise_time, fall_time, # Bazin bolometric
T_mid, T_amplitude, t_color, # sigmoid temperature: T (mid), T_amplitude, t_color
*baselines.values(), # per-band baselines
1.0, # reduced chi^2 placeholder
]
# Random observation times covering ±3 rise/fall windows
n_obs = 1000
t = np.sort(rng.uniform(reference_time - 3 * rise_time, reference_time + 3 * fall_time, n_obs))
band = rng.choice(list(band_wave_aa), size=n_obs)
# True model flux
flux = feature.model(t, band, *true_params)
# Poisson-like noise: S/N=10 at minimum flux
flux_err = np.sqrt(flux * np.min(flux)) / 10.0
flux = flux + rng.normal(0.0, flux_err)
print(f"Generated {n_obs} observations in bands {list(band_wave_aa)}")
Generated 1000 observations in bands ['g', 'r', 'i', 'z']
/home/runner/work/light-curve-python/light-curve-python/light-curve/light_curve/light_curve_py/warnings.py:19: ExperimentalWarning: Function light_curve.light_curve_py.features.rainbow.generic.RainbowFit is experimental and may cause any kind of troubles warn_experimental(message)
Fitting the light curve¶
Call feature(t, flux, sigma=flux_err, band=band). The result is a NumPy array of
fit parameters followed by the reduced χ² of the fit.
In [4]:
Copied!
result = feature(t, flux, sigma=flux_err, band=band)
# Parameter names are stored in feature.names
print("Parameter names:", feature.names)
print()
for name, val, true in zip(feature.names, result, true_params):
print(f" {name:25s} fitted={val:10.3f} true={true:10.3f}")
result = feature(t, flux, sigma=flux_err, band=band)
# Parameter names are stored in feature.names
print("Parameter names:", feature.names)
print()
for name, val, true in zip(feature.names, result, true_params):
print(f" {name:25s} fitted={val:10.3f} true={true:10.3f}")
Parameter names: ['reference_time', 'amplitude', 'rise_time', 'fall_time', 'T', 'T_amplitude', 't_color', 'baseline_g', 'baseline_r', 'baseline_i', 'baseline_z'] reference_time fitted= 59999.787 true= 60000.000 amplitude fitted= 0.968 true= 1.000 rise_time fitted= 5.332 true= 5.000 fall_time fitted= 31.953 true= 30.000 T fitted= 9554.275 true= 10000.000 T_amplitude fitted= 0.482 true= 0.500 t_color fitted= 10.303 true= 10.000 baseline_g fitted= 0.499 true= 0.504 baseline_r fitted= 0.588 true= 0.606 baseline_i fitted= 0.294 true= 0.306 baseline_z fitted= 0.286 true= 0.301
Visualising the fit¶
Plot the data and the best-fit model curve for each band.
In [5]:
Copied!
colors = {"g": "#2ca02c", "r": "#d62728", "i": "#9467bd", "z": "#8c564b"}
fig, ax = plt.subplots(figsize=(9, 4))
t_grid = np.linspace(t.min(), t.max(), 400)
for b in band_wave_aa:
mask = band == b
ax.errorbar(
t[mask], flux[mask], yerr=flux_err[mask],
fmt=".", color=colors[b], alpha=0.4, label=f"{b} data",
)
band_grid = np.full(len(t_grid), b)
ax.plot(t_grid, feature.model(t_grid, band_grid, *result), color=colors[b], lw=2)
ax.set_xlabel("Time (MJD)")
ax.set_ylabel("Flux density")
ax.set_title("RainbowFit — Bazin bolometric, sigmoid temperature")
ax.legend(ncol=2, fontsize=8)
plt.tight_layout()
plt.show()
colors = {"g": "#2ca02c", "r": "#d62728", "i": "#9467bd", "z": "#8c564b"}
fig, ax = plt.subplots(figsize=(9, 4))
t_grid = np.linspace(t.min(), t.max(), 400)
for b in band_wave_aa:
mask = band == b
ax.errorbar(
t[mask], flux[mask], yerr=flux_err[mask],
fmt=".", color=colors[b], alpha=0.4, label=f"{b} data",
)
band_grid = np.full(len(t_grid), b)
ax.plot(t_grid, feature.model(t_grid, band_grid, *result), color=colors[b], lw=2)
ax.set_xlabel("Time (MJD)")
ax.set_ylabel("Flux density")
ax.set_title("RainbowFit — Bazin bolometric, sigmoid temperature")
ax.legend(ncol=2, fontsize=8)
plt.tight_layout()
plt.show()
Parameter uncertainties¶
fit_and_get_errors returns a (params, errors) tuple. The errors are estimated by
iminuit from the Hessian of the chi-squared surface.
In [6]:
Copied!
params, errors = feature.fit_and_get_errors(t, flux, sigma=flux_err, band=band)
print(f"{'Parameter':25s} {'Value':>10s} {'Error':>10s} {'True':>10s}")
print("-" * 62)
for name, val, err, true in zip(feature.names, params, errors, true_params):
print(f"{name:25s} {val:10.3f} {err:10.3f} {true:10.3f}")
params, errors = feature.fit_and_get_errors(t, flux, sigma=flux_err, band=band)
print(f"{'Parameter':25s} {'Value':>10s} {'Error':>10s} {'True':>10s}")
print("-" * 62)
for name, val, err, true in zip(feature.names, params, errors, true_params):
print(f"{name:25s} {val:10.3f} {err:10.3f} {true:10.3f}")
Parameter Value Error True -------------------------------------------------------------- reference_time 59999.787 0.588 60000.000 amplitude 0.968 0.021 1.000 rise_time 5.332 0.232 5.000 fall_time 31.953 0.882 30.000 T 9554.275 481.343 10000.000 T_amplitude 0.482 0.020 0.500 t_color 10.303 0.951 10.000 baseline_g 0.499 0.036 0.504 baseline_r 0.588 0.037 0.606 baseline_i 0.294 0.034 0.306 baseline_z 0.286 0.036 0.301
Peak time¶
The bolometric peak time can differ from reference_time for asymmetric light curves.
feature.peak_time(params) computes it analytically from the bolometric term.
In [7]:
Copied!
t_peak = feature.peak_time(params)
print(f"Fitted peak time : {t_peak:.2f} MJD")
print(f"True peak time : {feature.peak_time(np.array(true_params)):.2f} MJD")
t_peak = feature.peak_time(params)
print(f"Fitted peak time : {t_peak:.2f} MJD")
print(f"True peak time : {feature.peak_time(np.array(true_params)):.2f} MJD")
Fitted peak time : 60007.97 MJD True peak time : 60007.68 MJD
Choosing different models¶
RainbowFit accepts pluggable bolometric, temperature, and spectral terms.
Bolometric terms¶
| Name | Parameters | Description |
|---|---|---|
'bazin' (default) |
reference_time, amplitude, rise_time, fall_time |
Bazin 2011 rise-and-fall |
'sigmoid' |
reference_time, amplitude, rise_time |
Logistic rise only |
'linexp' |
reference_time, amplitude, rise_time |
Linear rise, exponential decay |
'doublexp' |
reference_time, amplitude, time1, time2, p |
Exponential rise, decay with two exponentials |
Temperature terms¶
The sigmoid terms use the mid temperature T = $(T_\text{min}+T_\text{max})/2$ and a relative
amplitude T_amplitude = $(T_\text{max}-T_\text{min})/(T_\text{max}+T_\text{min})$
(T_amplitude = 0 is a constant temperature; a weak prior anchors it there unless the data demand a
temperature change). Strictly positive temperature corresponds to the independent box T > 0,
-1 < T_amplitude < 1.
| Name | Parameters | Description |
|---|---|---|
'sigmoid' (default) |
T, T_amplitude, t_color |
Temperature drops sigmoidally between $T(1+\texttt{T\_amplitude})$ and $T(1-\texttt{T\_amplitude})$ |
'constant' |
T |
Fixed temperature |
'delayed_sigmoid' |
T, T_amplitude, t_color, t_delay |
Sigmoid delayed relative to the bolometric peak |
Spectral terms¶
'planck' is a pure blackbody. The other terms describe SEDs that deviate from a blackbody
(UV blanketing or power-law spectra); each reduces to Planck at a null parameter value and
carries a weak Gaussian prior anchoring it there, so a true blackbody is recovered unbiased.
| Name | Parameters | Description |
|---|---|---|
'planck' (default) |
(none) | Standard blackbody $B_\nu(T)$ — no extra parameters |
'blanketed' |
lambda_scale |
UV-extincted blackbody $B_\nu(T)\,e^{-\tau}$, $\tau = I\,e^{-\lambda/\lambda_s}$; the extinction reach $\lambda_s$ is anchored to the characteristic temperature T (shared with the temperature term) so the blanketing depth stays fixed as the source cools |
'modified_bb' |
beta |
Planck tilted by a power law, $B_\nu(T)\,(\lambda/\lambda_\text{ref})^{\beta}$; $\beta = 0 \Rightarrow$ Planck |
'logparabola' |
sp_a, sp_b |
Planck × $e^{a L + b L^2}$, $L = \ln(\lambda/\lambda_\text{ref})$; $a = b = 0 \Rightarrow$ Planck |
'genwien' |
spec_k |
Generalized Wien $B_\nu \propto \nu^3 e^{-x^{\,\texttt{spec\_k}}}$, $x = h\nu/k_B T$; $\texttt{spec\_k} = 1 \Rightarrow$ Wien tail |
Both backends accept an optimizer argument: "iminuit" (default, robust Migrad) or
"least_squares" (scipy Trust Region Reflective, sharing the same analytic Jacobian — usually
faster, with an automatic fall-back to iminuit when needed).
In [8]:
Copied!
# Constant temperature — fewer parameters, useful for bluer transients
feature_const_T = RainbowFit.from_angstrom(
band_wave_aa,
with_baseline=False,
bolometric="bazin",
temperature="constant",
)
result_const = feature_const_T(t, flux, sigma=flux_err, band=band)
print("Parameters (constant T model):")
for name, val in zip(feature_const_T.names, result_const):
print(f" {name:25s} {val:.3f}")
# Constant temperature — fewer parameters, useful for bluer transients
feature_const_T = RainbowFit.from_angstrom(
band_wave_aa,
with_baseline=False,
bolometric="bazin",
temperature="constant",
)
result_const = feature_const_T(t, flux, sigma=flux_err, band=band)
print("Parameters (constant T model):")
for name, val in zip(feature_const_T.names, result_const):
print(f" {name:25s} {val:.3f}")
Parameters (constant T model): reference_time 59994.357 amplitude 2.007 rise_time 8.794 fall_time 87.390 T 9255.715
In [9]:
Copied!
# Blanketed blackbody — adds UV extinction parameters: lambda_scale
feature_blanketed = RainbowFit.from_angstrom(
band_wave_aa,
with_baseline=False,
bolometric="bazin",
temperature="sigmoid",
spectral="blanketed",
)
result_blanketed = feature_blanketed(t, flux, sigma=flux_err, band=band)
print("Parameters (blanketed spectral model):")
for name, val in zip(feature_blanketed.names, result_blanketed):
print(f" {name:25s} {val:.3f}")
# Blanketed blackbody — adds UV extinction parameters: lambda_scale
feature_blanketed = RainbowFit.from_angstrom(
band_wave_aa,
with_baseline=False,
bolometric="bazin",
temperature="sigmoid",
spectral="blanketed",
)
result_blanketed = feature_blanketed(t, flux, sigma=flux_err, band=band)
print("Parameters (blanketed spectral model):")
for name, val in zip(feature_blanketed.names, result_blanketed):
print(f" {name:25s} {val:.3f}")
Parameters (blanketed spectral model): reference_time 59981.092 T 20771.006 amplitude 3.671 rise_time 10.122 fall_time 85.846 T_amplitude 0.436 t_color 9.105 lambda_scale 0.293
Non-blackbody SED terms and the fast optimizer¶
Beyond 'blanketed', the spectral terms 'modified_bb', 'logparabola', and 'genwien'
describe SEDs that depart from a pure blackbody (a power-law tilt, a log-parabola, or a
generalized-Wien blue cutoff). Each reduces to Planck at a null parameter value
(beta = 0, sp_a = sp_b = 0, spec_k = 1), and a weak prior anchors it there — so fitting
them to genuinely blackbody data just recovers that null rather than inventing a deviation.
The fit itself can run on scipy's least_squares backend via optimizer="least_squares",
which shares the analytic Jacobian of the default iminuit path and is usually faster.
In [10]:
Copied!
# Each deviation term reduces to a pure blackbody at a "null" parameter value, with a weak
# prior anchoring it there. Fitting them to the (blackbody) data above recovers that null —
# i.e. no spurious deviation is introduced.
shape_params = {"modified_bb": ["beta"], "logparabola": ["sp_a", "sp_b"], "genwien": ["spec_k"]}
for spectral, shape in shape_params.items():
feat = RainbowFit.from_angstrom(
band_wave_aa, with_baseline=True,
bolometric="bazin", temperature="sigmoid", spectral=spectral,
)
res = feat(t, flux, sigma=flux_err, band=band)
shape_vals = ", ".join(f"{p}={res[feat.names.index(p)]:+.3f}" for p in shape)
print(f"{spectral:12s} reduced chi2 = {res[-1]:.2f} {shape_vals}")
# The same model fitted with the scipy least_squares backend (shares the analytic Jacobian).
feature_ls = RainbowFit.from_angstrom(
band_wave_aa, with_baseline=True,
bolometric="bazin", temperature="sigmoid", optimizer="least_squares",
)
result_ls = feature_ls(t, flux, sigma=flux_err, band=band)
print(f"\n{'least_squares':12s} reduced chi2 = {result_ls[-1]:.2f} (matches the default iminuit fit)")
# Each deviation term reduces to a pure blackbody at a "null" parameter value, with a weak
# prior anchoring it there. Fitting them to the (blackbody) data above recovers that null —
# i.e. no spurious deviation is introduced.
shape_params = {"modified_bb": ["beta"], "logparabola": ["sp_a", "sp_b"], "genwien": ["spec_k"]}
for spectral, shape in shape_params.items():
feat = RainbowFit.from_angstrom(
band_wave_aa, with_baseline=True,
bolometric="bazin", temperature="sigmoid", spectral=spectral,
)
res = feat(t, flux, sigma=flux_err, band=band)
shape_vals = ", ".join(f"{p}={res[feat.names.index(p)]:+.3f}" for p in shape)
print(f"{spectral:12s} reduced chi2 = {res[-1]:.2f} {shape_vals}")
# The same model fitted with the scipy least_squares backend (shares the analytic Jacobian).
feature_ls = RainbowFit.from_angstrom(
band_wave_aa, with_baseline=True,
bolometric="bazin", temperature="sigmoid", optimizer="least_squares",
)
result_ls = feature_ls(t, flux, sigma=flux_err, band=band)
print(f"\n{'least_squares':12s} reduced chi2 = {result_ls[-1]:.2f} (matches the default iminuit fit)")
modified_bb reduced chi2 = 0.50 beta=-0.125
logparabola reduced chi2 = 0.50 sp_a=-0.182, sp_b=+0.105
genwien reduced chi2 = 0.50 spec_k=+0.907
least_squares reduced chi2 = 0.50 (matches the default iminuit fit)
See also¶
- RainbowFit API reference
- Feature table — all extractors
- Parametric fits tutorial — single-band transient fits
In [ ]:
Copied!