{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Single-band parametric fits: BazinFit, LinexpFit, VillarFit\n",
    "\n",
    "[`BazinFit`](../api/fitting/#light_curve.BazinFit), [`LinexpFit`](../api/fitting/#light_curve.LinexpFit), and\n",
    "[`VillarFit`](../api/fitting/#light_curve.VillarFit) extract physically motivated parameters\n",
    "from transient **flux** light curves.\n",
    "All three operate on flux (not magnitude) — use positive flux values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %pip install light-curve"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "## Synthetic SN-like light curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import light_curve as licu\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "rng = np.random.default_rng(1)\n",
    "t = np.linspace(-20, 100, 80)\n",
    "\n",
    "# Ground truth: BazinFit functional form\n",
    "A, t0, rise, fall, B = 100.0, 10.0, 5.0, 25.0, 5.0\n",
    "true_params = np.array([A, B, t0, rise, fall])\n",
    "\n",
    "bazin_ref = licu.BazinFit('lmsder')  # used only for .model()\n",
    "flux_true = bazin_ref.model(t, true_params)\n",
    "flux_err  = np.sqrt(np.abs(flux_true)) * 0.1\n",
    "flux      = flux_true + rng.normal(0, flux_err)\n",
    "print(f'n_obs: {len(t)},  peak flux: {flux.max():.1f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "## BazinFit\n",
    "\n",
    "[`BazinFit`](../api/fitting/#light_curve.BazinFit) fits a 5-parameter rising/falling exponential (Bazin et al. 2009):\n",
    "\n",
    "$$f(t) = A \\cdot \\frac{e^{-(t-t_0)/t_\\text{fall}}}{1 + e^{-(t-t_0)/t_\\text{rise}}} + B$$\n",
    "\n",
    "Output: amplitude $A$, baseline $B$, reference time $t_0$, rise time $\\tau_\\text{rise}$,\n",
    "fall time $\\tau_\\text{fall}$, and reduced $\\chi^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import light_curve as licu\n",
    "import numpy as np\n",
    "\n",
    "rng = np.random.default_rng(1)\n",
    "t = np.linspace(-20, 100, 80)\n",
    "A, t0, rise, fall, B = 100.0, 10.0, 5.0, 25.0, 5.0\n",
    "true_params = np.array([A, B, t0, rise, fall])\n",
    "bazin_ref = licu.BazinFit('lmsder')\n",
    "flux_true = bazin_ref.model(t, true_params)\n",
    "flux_err  = np.sqrt(np.abs(flux_true)) * 0.1\n",
    "flux      = flux_true + rng.normal(0, flux_err)\n",
    "\n",
    "bazin = licu.BazinFit(algorithm='mcmc-ceres')\n",
    "result_b = bazin(t, flux, flux_err)\n",
    "for name, val in zip(bazin.names, result_b):\n",
    "    print(f'  {name:35s} = {val:.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "## LinexpFit\n",
    "\n",
    "[`LinexpFit`](../api/fitting/#light_curve.LinexpFit) fits a linear-times-exponential (Karpenka et al. 2012) — simple rising part, fast:\n",
    "\n",
    "$$f(t) = A \\cdot (1 + B\\,(t - t_0)) \\cdot e^{-(t - t_0) / t_\\text{fall}} \\cdot \\mathbf{1}_{t > t_0} + C$$\n",
    "\n",
    "Output: amplitude, slope $B$, reference time, fall time, baseline, reduced $\\chi^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import light_curve as licu\n",
    "import numpy as np\n",
    "\n",
    "rng = np.random.default_rng(1)\n",
    "t = np.linspace(-20, 100, 80)\n",
    "A, t0, rise, fall, B = 100.0, 10.0, 5.0, 25.0, 5.0\n",
    "bazin_ref = licu.BazinFit('lmsder')\n",
    "flux_true = bazin_ref.model(t, np.array([A, B, t0, rise, fall]))\n",
    "flux_err  = np.sqrt(np.abs(flux_true)) * 0.1\n",
    "flux      = flux_true + rng.normal(0, flux_err)\n",
    "\n",
    "linexp = licu.LinexpFit(algorithm='ceres')\n",
    "result_l = linexp(t, flux, flux_err)\n",
    "for name, val in zip(linexp.names, result_l):\n",
    "    print(f'  {name:35s} = {val:.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "## VillarFit\n",
    "\n",
    "[`VillarFit`](../api/fitting/#light_curve.VillarFit) fits the Villar et al. 2019 function — Gaussian-plateau model for SN classification:\n",
    "\n",
    "$$f(t) = \\frac{A + \\beta\\,(t - t_0)}{1 + e^{-(t - t_0)/t_\\text{rise}}} \\cdot\n",
    "\\begin{cases}\n",
    "1 & t \\leq t_0 + t_\\text{plateau} \\\\\n",
    "e^{-(t - t_0 - t_\\text{plateau})/t_\\text{fall}} & t > t_0 + t_\\text{plateau}\n",
    "\\end{cases} + B$$\n",
    "\n",
    "Output: amplitude $A$, baseline $B$, reference time $t_0$, rise time, fall time,\n",
    "plateau slope $\\beta$, plateau duration, reduced $\\chi^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import light_curve as licu\n",
    "import numpy as np\n",
    "\n",
    "rng = np.random.default_rng(1)\n",
    "t = np.linspace(-20, 100, 80)\n",
    "A, t0, rise, fall, B = 100.0, 10.0, 5.0, 25.0, 5.0\n",
    "bazin_ref = licu.BazinFit('lmsder')\n",
    "flux_true = bazin_ref.model(t, np.array([A, B, t0, rise, fall]))\n",
    "flux_err  = np.sqrt(np.abs(flux_true)) * 0.1\n",
    "flux      = flux_true + rng.normal(0, flux_err)\n",
    "\n",
    "villar = licu.VillarFit(algorithm='ceres')\n",
    "result_v = villar(t, flux, flux_err)\n",
    "for name, val in zip(villar.names, result_v):\n",
    "    print(f'  {name:35s} = {val:.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "## Comparing the fits\n",
    "\n",
    "`.model(t, params)` evaluates the fitted curve on any time grid — no need to\n",
    "re-implement the equations. Pass the full parameter array; the method uses only\n",
    "the first N entries (ignoring the trailing reduced χ²)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "import light_curve as licu\n",
    "import numpy as np\n",
    "\n",
    "rng = np.random.default_rng(1)\n",
    "t = np.linspace(-20, 100, 80)\n",
    "A, t0, rise, fall, B = 100.0, 10.0, 5.0, 25.0, 5.0\n",
    "bazin_ref = licu.BazinFit('lmsder')\n",
    "flux_true = bazin_ref.model(t, np.array([A, B, t0, rise, fall]))\n",
    "flux_err  = np.sqrt(np.abs(flux_true)) * 0.1\n",
    "flux      = flux_true + rng.normal(0, flux_err)\n",
    "\n",
    "bazin  = licu.BazinFit('mcmc-ceres')\n",
    "linexp = licu.LinexpFit('ceres')\n",
    "villar = licu.VillarFit('ceres')\n",
    "\n",
    "result_b = bazin(t,  flux, flux_err)\n",
    "result_l = linexp(t, flux, flux_err)\n",
    "result_v = villar(t, flux, flux_err)\n",
    "\n",
    "t_grid = np.linspace(t.min(), t.max(), 300)\n",
    "flux_b = bazin.model(t_grid,  result_b)\n",
    "flux_l = linexp.model(t_grid, result_l)\n",
    "flux_v = villar.model(t_grid, result_v)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9, 4))\n",
    "ax.errorbar(t, flux, yerr=flux_err, fmt='.', color='gray', alpha=0.5, label='data')\n",
    "ax.plot(t_grid, flux_b, label=f'BazinFit  χ²={result_b[-1]:.2f}', lw=2)\n",
    "ax.plot(t_grid, flux_l, label=f'LinexpFit χ²={result_l[-1]:.2f}', lw=2, ls='--')\n",
    "ax.plot(t_grid, flux_v, label=f'VillarFit χ²={result_v[-1]:.2f}', lw=2, ls=':')\n",
    "ax.set_xlabel('Time (days)')\n",
    "ax.set_ylabel('Flux')\n",
    "ax.legend()\n",
    "ax.set_title('Parametric fit comparison')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "metadata": {},
   "source": "## Solvers\n\nAll three classes accept an `algorithm` argument. All return the same 6 outputs\n(5 fit parameters + reduced χ²) regardless of solver.\n\nSingle-stage solvers go straight to a local minimum:\n\n| Algorithm | Speed | Notes |\n|-----------|-------|-------|\n| `'ceres'` | Fast | Gradient-based via Google Ceres — **preferred** |\n| `'lmsder'` | Fast | Levenberg–Marquardt via GSL; not available on Windows |\n| `'nuts'` | Slow | No-U-Turn Sampler MCMC — **preferred** |\n| `'mcmc'` | Slow | Random-walk MCMC |\n\nTwo-stage solvers (dash-separated) **chain a sampler into a gradient descent**:\nthe sampler (`nuts` or `mcmc`) runs first to broadly explore parameter space\nand find the region of the global minimum, then the gradient-based refiner\n(`ceres` or `lmsder`) starts from that best point and converges precisely.\nThis combines global search robustness with fast local refinement:\n\n| Algorithm | Stage 1 | Stage 2 | Notes |\n|-----------|---------|---------|-------|\n| `'nuts-ceres'` | NUTS | Ceres | **Recommended**: most robust |\n| `'mcmc-ceres'` | MCMC | Ceres | Faster stage 1, slightly less thorough |\n| `'nuts-lmsder'` | NUTS | LM | Not available on Windows |\n| `'mcmc-lmsder'` | MCMC | LM | Not available on Windows |\n\nUse a single-stage `'ceres'` for speed-critical batch processing when light curves\nare well-sampled. Use `'nuts-ceres'` when data are sparse or noisy and gradient\ndescent is likely to get stuck.\n\n**Solver settings** (keyword arguments to the constructor):\n\n| Parameter | Default | Applies to | Effect |\n|-----------|---------|------------|--------|\n| `mcmc_niter` | 128 | `mcmc*` | MCMC chain length — increase for better coverage |\n| `ceres_niter` | 10 | `*ceres` | Ceres refinement iterations |\n| `lmsder_niter` | 10 | `*lmsder` | LM refinement iterations |\n| `ceres_loss_reg` | `None` | `*ceres` | Huber loss threshold for outlier rejection |\n| `init` | `None` | `mcmc*` | Initial parameter guess (list of 5 floats or `None`s) |\n| `bounds` | `None` | `mcmc*` | Parameter bounds (list of `(lo, hi)` pairs or `None`s) |"
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": "## See also\n\n- All three features work on **flux** light curves only (not magnitude).\n- For linear trend and intercept, see [`LinearFit`](../api/linear/#light_curve.LinearFit) and [`LinearTrend`](../api/linear/#light_curve.LinearTrend).\n- For multi-band transient fitting with a physical blackbody model, see the\n  [Rainbow tutorial](../rainbow/).\n- [API reference](../api/fitting/)"
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}