Single-band parametric fits: BazinFit, LinexpFit, VillarFit¶

BazinFit, LinexpFit, and VillarFit extract physically motivated parameters from transient flux light curves. All three operate on flux (not magnitude) — use positive flux values.

In [1]:

Copied!

# %pip install light-curve
# %pip install light-curve

Synthetic SN-like light curve¶

In [2]:

Copied!





import light_curve as licu
import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(1)
t = np.linspace(-20, 100, 80)

# Ground truth: BazinFit functional form
A, t0, rise, fall, B = 100.0, 10.0, 5.0, 25.0, 5.0
true_params = np.array([A, B, t0, rise, fall])

bazin_ref = licu.BazinFit('lmsder')  # used only for .model()
flux_true = bazin_ref.model(t, true_params)
flux_err  = np.sqrt(np.abs(flux_true)) * 0.1
flux      = flux_true + rng.normal(0, flux_err)
print(f'n_obs: {len(t)},  peak flux: {flux.max():.1f}')
import light_curve as licu
import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(1)
t = np.linspace(-20, 100, 80)

# Ground truth: BazinFit functional form
A, t0, rise, fall, B = 100.0, 10.0, 5.0, 25.0, 5.0
true_params = np.array([A, B, t0, rise, fall])

bazin_ref = licu.BazinFit('lmsder')  # used only for .model()
flux_true = bazin_ref.model(t, true_params)
flux_err  = np.sqrt(np.abs(flux_true)) * 0.1
flux      = flux_true + rng.normal(0, flux_err)
print(f'n_obs: {len(t)},  peak flux: {flux.max():.1f}')

n_obs: 80,  peak flux: 65.6

BazinFit¶

BazinFit fits a 5-parameter rising/falling exponential (Bazin et al. 2009):

$$f(t) = A \cdot \frac{e^{-(t-t_0)/t_\text{fall}}}{1 + e^{-(t-t_0)/t_\text{rise}}} + B$$

Output: amplitude $A$, baseline $B$, reference time $t_0$, rise time $\tau_\text{rise}$, fall time $\tau_\text{fall}$, and reduced $\chi^2$.

In [3]:

Copied!





import light_curve as licu
import numpy as np

rng = np.random.default_rng(1)
t = np.linspace(-20, 100, 80)
A, t0, rise, fall, B = 100.0, 10.0, 5.0, 25.0, 5.0
true_params = np.array([A, B, t0, rise, fall])
bazin_ref = licu.BazinFit('lmsder')
flux_true = bazin_ref.model(t, true_params)
flux_err  = np.sqrt(np.abs(flux_true)) * 0.1
flux      = flux_true + rng.normal(0, flux_err)

bazin = licu.BazinFit(algorithm='mcmc-ceres')
result_b = bazin(t, flux, flux_err)
for name, val in zip(bazin.names, result_b):
    print(f'  {name:35s} = {val:.4f}')
import light_curve as licu
import numpy as np

rng = np.random.default_rng(1)
t = np.linspace(-20, 100, 80)
A, t0, rise, fall, B = 100.0, 10.0, 5.0, 25.0, 5.0
true_params = np.array([A, B, t0, rise, fall])
bazin_ref = licu.BazinFit('lmsder')
flux_true = bazin_ref.model(t, true_params)
flux_err  = np.sqrt(np.abs(flux_true)) * 0.1
flux      = flux_true + rng.normal(0, flux_err)

bazin = licu.BazinFit(algorithm='mcmc-ceres')
result_b = bazin(t, flux, flux_err)
for name, val in zip(bazin.names, result_b):
    print(f'  {name:35s} = {val:.4f}')

  bazin_fit_amplitude                 = 100.0255
  bazin_fit_baseline                  = 5.0111
  bazin_fit_reference_time            = 10.0869
  bazin_fit_rise_time                 = 5.0271
  bazin_fit_fall_time                 = 24.8894
  bazin_fit_reduced_chi2              = 0.7394

LinexpFit¶

LinexpFit fits a linear-times-exponential (Karpenka et al. 2012) — simple rising part, fast:

$$f(t) = A \cdot (1 + B\,(t - t_0)) \cdot e^{-(t - t_0) / t_\text{fall}} \cdot \mathbf{1}_{t > t_0} + C$$

Output: amplitude, slope $B$, reference time, fall time, baseline, reduced $\chi^2$.

In [4]:

Copied!





import light_curve as licu
import numpy as np

rng = np.random.default_rng(1)
t = np.linspace(-20, 100, 80)
A, t0, rise, fall, B = 100.0, 10.0, 5.0, 25.0, 5.0
bazin_ref = licu.BazinFit('lmsder')
flux_true = bazin_ref.model(t, np.array([A, B, t0, rise, fall]))
flux_err  = np.sqrt(np.abs(flux_true)) * 0.1
flux      = flux_true + rng.normal(0, flux_err)

linexp = licu.LinexpFit(algorithm='ceres')
result_l = linexp(t, flux, flux_err)
for name, val in zip(linexp.names, result_l):
    print(f'  {name:35s} = {val:.4f}')
import light_curve as licu
import numpy as np

rng = np.random.default_rng(1)
t = np.linspace(-20, 100, 80)
A, t0, rise, fall, B = 100.0, 10.0, 5.0, 25.0, 5.0
bazin_ref = licu.BazinFit('lmsder')
flux_true = bazin_ref.model(t, np.array([A, B, t0, rise, fall]))
flux_err  = np.sqrt(np.abs(flux_true)) * 0.1
flux      = flux_true + rng.normal(0, flux_err)

linexp = licu.LinexpFit(algorithm='ceres')
result_l = linexp(t, flux, flux_err)
for name, val in zip(linexp.names, result_l):
    print(f'  {name:35s} = {val:.4f}')

  linexp_fit_amplitude                = 199.1246
  linexp_fit_reference_time           = -29.9739
  linexp_fit_fall_time                = 50.4585
  linexp_fit_baseline                 = -36.0082
  linexp_fit_reduced_chi2             = 333.0247

VillarFit¶

VillarFit fits the Villar et al. 2019 function — Gaussian-plateau model for SN classification:

$$f(t) = \frac{A + \beta\,(t - t_0)}{1 + e^{-(t - t_0)/t_\text{rise}}} \cdot \begin{cases} 1 & t \leq t_0 + t_\text{plateau} \\ e^{-(t - t_0 - t_\text{plateau})/t_\text{fall}} & t > t_0 + t_\text{plateau} \end{cases} + B$$

Output: amplitude $A$, baseline $B$, reference time $t_0$, rise time, fall time, plateau slope $\beta$, plateau duration, reduced $\chi^2$.

In [5]:

Copied!





import light_curve as licu
import numpy as np

rng = np.random.default_rng(1)
t = np.linspace(-20, 100, 80)
A, t0, rise, fall, B = 100.0, 10.0, 5.0, 25.0, 5.0
bazin_ref = licu.BazinFit('lmsder')
flux_true = bazin_ref.model(t, np.array([A, B, t0, rise, fall]))
flux_err  = np.sqrt(np.abs(flux_true)) * 0.1
flux      = flux_true + rng.normal(0, flux_err)

villar = licu.VillarFit(algorithm='ceres')
result_v = villar(t, flux, flux_err)
for name, val in zip(villar.names, result_v):
    print(f'  {name:35s} = {val:.4f}')
import light_curve as licu
import numpy as np

rng = np.random.default_rng(1)
t = np.linspace(-20, 100, 80)
A, t0, rise, fall, B = 100.0, 10.0, 5.0, 25.0, 5.0
bazin_ref = licu.BazinFit('lmsder')
flux_true = bazin_ref.model(t, np.array([A, B, t0, rise, fall]))
flux_err  = np.sqrt(np.abs(flux_true)) * 0.1
flux      = flux_true + rng.normal(0, flux_err)

villar = licu.VillarFit(algorithm='ceres')
result_v = villar(t, flux, flux_err)
for name, val in zip(villar.names, result_v):
    print(f'  {name:35s} = {val:.4f}')

  villar_fit_amplitude                = 54.0211
  villar_fit_baseline                 = 5.8620
  villar_fit_reference_time           = 3.1129
  villar_fit_rise_time                = 4.6050
  villar_fit_fall_time                = 24.3447
  villar_fit_plateau_rel_amplitude    = 0.0000
  villar_fit_plateau_duration         = 18.1314
  villar_fit_reduced_chi2             = 25.9662

WARNING: Logging before InitGoogleLogging() is written to STDERR
E20260626 19:46:53.791675 139795735391104 program_evaluator.h:286] Accumulated cost = inf is not a finite number. Evaluation failed.
E20260626 19:46:53.792982 139795735391104 program_evaluator.h:286] Accumulated cost = inf is not a finite number. Evaluation failed.

Comparing the fits¶

.model(t, params) evaluates the fitted curve on any time grid — no need to re-implement the equations. Pass the full parameter array; the method uses only the first N entries (ignoring the trailing reduced χ²).

In [6]:

Copied!





import light_curve as licu
import numpy as np

rng = np.random.default_rng(1)
t = np.linspace(-20, 100, 80)
A, t0, rise, fall, B = 100.0, 10.0, 5.0, 25.0, 5.0
bazin_ref = licu.BazinFit('lmsder')
flux_true = bazin_ref.model(t, np.array([A, B, t0, rise, fall]))
flux_err  = np.sqrt(np.abs(flux_true)) * 0.1
flux      = flux_true + rng.normal(0, flux_err)

bazin  = licu.BazinFit('mcmc-ceres')
linexp = licu.LinexpFit('ceres')
villar = licu.VillarFit('ceres')

result_b = bazin(t,  flux, flux_err)
result_l = linexp(t, flux, flux_err)
result_v = villar(t, flux, flux_err)

t_grid = np.linspace(t.min(), t.max(), 300)
flux_b = bazin.model(t_grid,  result_b)
flux_l = linexp.model(t_grid, result_l)
flux_v = villar.model(t_grid, result_v)

fig, ax = plt.subplots(figsize=(9, 4))
ax.errorbar(t, flux, yerr=flux_err, fmt='.', color='gray', alpha=0.5, label='data')
ax.plot(t_grid, flux_b, label=f'BazinFit  χ²={result_b[-1]:.2f}', lw=2)
ax.plot(t_grid, flux_l, label=f'LinexpFit χ²={result_l[-1]:.2f}', lw=2, ls='--')
ax.plot(t_grid, flux_v, label=f'VillarFit χ²={result_v[-1]:.2f}', lw=2, ls=':')
ax.set_xlabel('Time (days)')
ax.set_ylabel('Flux')
ax.legend()
ax.set_title('Parametric fit comparison')
plt.tight_layout()
plt.show()
import light_curve as licu
import numpy as np

rng = np.random.default_rng(1)
t = np.linspace(-20, 100, 80)
A, t0, rise, fall, B = 100.0, 10.0, 5.0, 25.0, 5.0
bazin_ref = licu.BazinFit('lmsder')
flux_true = bazin_ref.model(t, np.array([A, B, t0, rise, fall]))
flux_err  = np.sqrt(np.abs(flux_true)) * 0.1
flux      = flux_true + rng.normal(0, flux_err)

bazin  = licu.BazinFit('mcmc-ceres')
linexp = licu.LinexpFit('ceres')
villar = licu.VillarFit('ceres')

result_b = bazin(t,  flux, flux_err)
result_l = linexp(t, flux, flux_err)
result_v = villar(t, flux, flux_err)

t_grid = np.linspace(t.min(), t.max(), 300)
flux_b = bazin.model(t_grid,  result_b)
flux_l = linexp.model(t_grid, result_l)
flux_v = villar.model(t_grid, result_v)

fig, ax = plt.subplots(figsize=(9, 4))
ax.errorbar(t, flux, yerr=flux_err, fmt='.', color='gray', alpha=0.5, label='data')
ax.plot(t_grid, flux_b, label=f'BazinFit  χ²={result_b[-1]:.2f}', lw=2)
ax.plot(t_grid, flux_l, label=f'LinexpFit χ²={result_l[-1]:.2f}', lw=2, ls='--')
ax.plot(t_grid, flux_v, label=f'VillarFit χ²={result_v[-1]:.2f}', lw=2, ls=':')
ax.set_xlabel('Time (days)')
ax.set_ylabel('Flux')
ax.legend()
ax.set_title('Parametric fit comparison')
plt.tight_layout()
plt.show()

E20260626 19:46:53.842248 139795735391104 program_evaluator.h:286] Accumulated cost = inf is not a finite number. Evaluation failed.
E20260626 19:46:53.843553 139795735391104 program_evaluator.h:286] Accumulated cost = inf is not a finite number. Evaluation failed.

No description has been provided for this image

Solvers¶

All three classes accept an algorithm argument. All return the same 6 outputs (5 fit parameters + reduced χ²) regardless of solver.

Single-stage solvers go straight to a local minimum:

Algorithm	Speed	Notes
`'ceres'`	Fast	Gradient-based via Google Ceres — preferred
`'lmsder'`	Fast	Levenberg–Marquardt via GSL; not available on Windows
`'nuts'`	Slow	No-U-Turn Sampler MCMC — preferred
`'mcmc'`	Slow	Random-walk MCMC

Two-stage solvers (dash-separated) chain a sampler into a gradient descent: the sampler (nuts or mcmc) runs first to broadly explore parameter space and find the region of the global minimum, then the gradient-based refiner (ceres or lmsder) starts from that best point and converges precisely. This combines global search robustness with fast local refinement:

Algorithm	Stage 1	Stage 2	Notes
`'nuts-ceres'`	NUTS	Ceres	Recommended: most robust
`'mcmc-ceres'`	MCMC	Ceres	Faster stage 1, slightly less thorough
`'nuts-lmsder'`	NUTS	LM	Not available on Windows
`'mcmc-lmsder'`	MCMC	LM	Not available on Windows

Use a single-stage 'ceres' for speed-critical batch processing when light curves are well-sampled. Use 'nuts-ceres' when data are sparse or noisy and gradient descent is likely to get stuck.

Solver settings (keyword arguments to the constructor):

Parameter	Default	Applies to	Effect
`mcmc_niter`	128	`mcmc*`	MCMC chain length — increase for better coverage
`ceres_niter`	10	`*ceres`	Ceres refinement iterations
`lmsder_niter`	10	`*lmsder`	LM refinement iterations
`ceres_loss_reg`	`None`	`*ceres`	Huber loss threshold for outlier rejection
`init`	`None`	`mcmc*`	Initial parameter guess (list of 5 floats or `None`s)
`bounds`	`None`	`mcmc*`	Parameter bounds (list of `(lo, hi)` pairs or `None`s)