{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Periodogram\n",
    "\n",
    "The [`Periodogram`](../api/periodogram/#light_curve.Periodogram) feature extractor finds peaks in the **Lomb–Scargle\n",
    "periodogram** of an unevenly sampled time series, returning period and signal-to-noise\n",
    "ratio for each peak:\n",
    "\n",
    "$$\\mathrm{S/N}(\\omega_\\mathrm{peak}) = \\frac{P(\\omega_\\mathrm{peak}) - \\langle P(\\omega) \\rangle}{\\sigma_{P(\\omega)}}$$\n",
    "\n",
    "This tutorial covers parameter selection, reading the full power spectrum, multi-period\n",
    "detection, and using spectrum/phase features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %pip install light-curve matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "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(0)\n",
    "PERIOD = 23.5  # days\n",
    "t = np.sort(rng.uniform(0, 300, 150))\n",
    "m = 15.0 + 0.4 * np.sin(2 * np.pi * t / PERIOD) + rng.normal(0, 0.04, 150)\n",
    "print(f'n_obs: {len(t)}, time span: {t[-1] - t[0]:.0f} d, true period: {PERIOD} d')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "## Basic usage\n",
    "\n",
    "`Periodogram(peaks=N)` returns the `N` strongest periods and their S/N ratios.\n",
    "With `peaks=1` you get two numbers: `[period, snr]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "rng = np.random.default_rng(0)\n",
    "PERIOD = 23.5\n",
    "t = np.sort(rng.uniform(0, 300, 150))\n",
    "m = 15.0 + 0.4 * np.sin(2 * np.pi * t / PERIOD) + rng.normal(0, 0.04, 150)\n",
    "\n",
    "pg = licu.Periodogram(peaks=1, fast=True)\n",
    "result = pg(t, m)\n",
    "for name, val in zip(pg.names, result):\n",
    "    print(f'  {name}: {val:.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "## Key parameters\n",
    "\n",
    "| Parameter | Default | What it controls |\n",
    "|-----------|---------|------------------|\n",
    "| `peaks` | 1 | Number of period/S/N pairs returned |\n",
    "| `resolution` | 10 | Frequency grid oversampling factor |\n",
    "| `max_freq_factor` | 1.0 | Upper limit multiplier on the Nyquist frequency |\n",
    "| `nyquist` | `'average'` | How the Nyquist frequency is estimated (`'average'`, `'median'`, or a float quantile) |\n",
    "| `fast` | `True` | NFFT-based O(N log N) algorithm; `False` uses exact O(N²) |\n",
    "| `normalization` | `'psd'` | Power normalisation (see below) |\n",
    "| `freqs` | `None` | Explicit angular frequency grid; overrides `resolution`, `max_freq_factor`, `nyquist` |\n",
    "\n",
    "**`resolution`** controls how densely the frequency grid is sampled relative to the\n",
    "minimum resolvable frequency (1/T_baseline). Higher values find weaker signals but\n",
    "increase runtime linearly. Values of 10–20 work well for most survey data.\n",
    "\n",
    "**`normalization`** options:\n",
    "\n",
    "| Value | Range | Notes |\n",
    "|-------|-------|-------|\n",
    "| `'psd'` (default) | [0, ∞) | Raw power; good for S/N comparisons |\n",
    "| `'standard'` | [0, 1] | $P_\\text{std} = P \\cdot 2/(n-1)$; matches astropy |\n",
    "| `'model'` | [0, ∞) | $P_\\text{std}/(1-P_\\text{std})$; matches astropy |\n",
    "| `'log'` | [0, ∞) | $-\\ln(1-P_\\text{std})$; matches astropy |\n",
    "\n",
    "**`freqs`**: pass an explicit angular frequency array (radians per time unit) to\n",
    "control the grid exactly. For `fast=True`, the grid must be\n",
    "`np.linspace(0.0, max_freq, 2**k + 1)` for integer k."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "## Power spectrum\n",
    "\n",
    "[`Periodogram`](../api/periodogram/#light_curve.Periodogram)`.freq_power(t, m)` returns the full `(angular_frequencies, power)` arrays,\n",
    "so you can visualise or post-process the complete spectrum without any external library.\n",
    "Frequencies are in **radians per time unit**; convert to period with $T = 2\\pi/\\omega$."
   ]
  },
  {
   "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(0)\n",
    "PERIOD = 23.5\n",
    "t = np.sort(rng.uniform(0, 300, 150))\n",
    "m = 15.0 + 0.4 * np.sin(2 * np.pi * t / PERIOD) + rng.normal(0, 0.04, 150)\n",
    "\n",
    "pg = licu.Periodogram(peaks=1, resolution=10)\n",
    "omega, power = pg.freq_power(t, m)  # angular freqs (rad/day)\n",
    "\n",
    "# Skip omega=0 to avoid division by zero when converting to period\n",
    "mask = omega > 0\n",
    "periods = 2 * np.pi / omega[mask]\n",
    "pwr = power[mask]\n",
    "\n",
    "best_period = periods[np.argmax(pwr)]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9, 3))\n",
    "ax.plot(periods, pwr, lw=0.7, color='steelblue')\n",
    "ax.axvline(best_period, color='red', lw=1.5, ls='--', label=f'Best = {best_period:.2f} d')\n",
    "ax.axvline(PERIOD, color='green', lw=1.5, ls=':', label=f'True = {PERIOD} d')\n",
    "ax.set_xscale('log')\n",
    "ax.set_xlabel('Period (days)')\n",
    "ax.set_ylabel('LS power (psd)')\n",
    "ax.set_title('Lomb–Scargle power spectrum')\n",
    "ax.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "print(f'Recovered period: {best_period:.3f} d  (true: {PERIOD} d)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "## Multiple peaks and multi-period detection\n",
    "\n",
    "When the light curve contains several periodic signals, `peaks=N` returns the N\n",
    "strongest ones. S/N drops sharply for noise peaks, making it a natural threshold.\n",
    "Here we synthesise a light curve with **two real periods** and ask for five peaks:"
   ]
  },
  {
   "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(0)\n",
    "P1, P2 = 23.5, 7.3  # two embedded periods\n",
    "t = np.sort(rng.uniform(0, 400, 200))\n",
    "m = (15.0\n",
    "     + 0.40 * np.sin(2 * np.pi * t / P1)\n",
    "     + 0.25 * np.sin(2 * np.pi * t / P2)\n",
    "     + rng.normal(0, 0.04, 200))\n",
    "\n",
    "pg5 = licu.Periodogram(peaks=5, resolution=20)\n",
    "result5 = pg5(t, m)\n",
    "\n",
    "print(f'True periods: {P1} d and {P2} d')\n",
    "print('Top 5 peaks:')\n",
    "for i in range(5):\n",
    "    period = result5[i * 2]\n",
    "    snr = result5[i * 2 + 1]\n",
    "    marker = ' ← real' if abs(period - P1) < 1 or abs(period - P2) < 1 else ''\n",
    "    print(f'  peak {i + 1}: {period:6.2f} d, S/N = {snr:5.2f}{marker}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": "## Spectrum features\n\n`features=` applies feature extractors to the **power spectrum** treated as a time\nseries (frequency as time, power as magnitude). This captures the shape of the\nspectrum itself — useful for distinguishing sharp peaks from broad or noisy ones.\nHere we add [`Amplitude`](../api/variability/#light_curve.Amplitude), [`Kurtosis`](../api/variability/#light_curve.Kurtosis),\n[`Skew`](../api/variability/#light_curve.Skew), and [`BeyondNStd`](../api/variability/#light_curve.BeyondNStd):"
  },
  {
   "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(0)\n",
    "PERIOD = 23.5\n",
    "t = np.sort(rng.uniform(0, 300, 150))\n",
    "m = 15.0 + 0.4 * np.sin(2 * np.pi * t / PERIOD) + rng.normal(0, 0.04, 150)\n",
    "\n",
    "pg_spec = licu.Periodogram(\n",
    "    peaks=1,\n",
    "    features=[licu.Amplitude(), licu.Kurtosis(), licu.Skew(), licu.BeyondNStd(nstd=2)],\n",
    ")\n",
    "res_spec = pg_spec(t, m)\n",
    "print('Periodogram + spectrum features:')\n",
    "for name, val in zip(pg_spec.names, res_spec):\n",
    "    print(f'  {name:50s} = {val:.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "metadata": {},
   "source": [
    "## Phase-folded features\n",
    "\n",
    "`phase_features=` applies extractors to the light curve **folded at the best period**.\n",
    "The phase runs 0 → 1 (phase 0 = magnitude minimum). Features designed for\n",
    "smoothly varying series react very differently to a folded periodic signal versus\n",
    "the original unordered time series.\n",
    "\n",
    "The **[`LaflerKinmanStringLength`](../api/variability/#light_curve.LaflerKinmanStringLength)** $\\theta$ is the normalised sum of squared\n",
    "successive differences of phase-sorted magnitudes — a clean, smooth periodic curve\n",
    "gives $\\theta \\ll 1$, while a noisy or aperiodic series gives $\\theta \\approx 1$.\n",
    "Similarly, **[`EtaE`](../api/variability/#light_curve.EtaE)** (von Neumann $\\eta^e$) is small for a random series and large\n",
    "for a phase-sorted periodic one, making it sensitive to regularity.\n",
    "\n",
    "Extracting these *without* phase folding returns uninformative values; the table\n",
    "below shows the dramatic difference:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "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(0)\n",
    "PERIOD = 23.5\n",
    "t = np.sort(rng.uniform(0, 300, 150))\n",
    "m = 15.0 + 0.4 * np.sin(2 * np.pi * t / PERIOD) + rng.normal(0, 0.04, 150)\n",
    "\n",
    "pg_phase = licu.Periodogram(\n",
    "    peaks=1,\n",
    "    phase_features=[licu.LaflerKinmanStringLength(), licu.EtaE()],\n",
    ")\n",
    "result_pf = pg_phase(t, m)\n",
    "best_p = result_pf[0]\n",
    "\n",
    "# Compare with unfolded values\n",
    "lk = licu.LaflerKinmanStringLength()\n",
    "etae = licu.EtaE()\n",
    "print(f'{'Feature':<35} {'non-folded':>12} {'phase-folded':>14}')\n",
    "print('-' * 63)\n",
    "print(f'{'string_length':<35} {lk(t, m)[0]:>12.4f} {result_pf[2]:>14.4f}')\n",
    "print(f'{'eta_e':<35} {etae(t, m)[0]:>12.4f} {result_pf[3]:>14.4f}')\n",
    "print()\n",
    "print(f'Recovered period: {best_p:.3f} d  (true: {PERIOD} d)')\n",
    "\n",
    "# Plot the phase-folded light curve\n",
    "phase = (t % best_p) / best_p\n",
    "fig, ax = plt.subplots(figsize=(6, 3))\n",
    "ax.plot(phase, m, '.', alpha=0.5)\n",
    "ax.set_xlabel('Phase')\n",
    "ax.set_ylabel('Magnitude')\n",
    "ax.invert_yaxis()\n",
    "plt.xlim(0, 1)\n",
    "ax.set_title(f'Phase-folded light curve  (P = {best_p:.2f} d)')\n",
    "plt.tight_layout()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14",
   "metadata": {},
   "source": [
    "## Explicit frequency grid\n",
    "\n",
    "Pass `freqs=` to fix the angular frequency grid (radians per time unit) instead of\n",
    "computing it from the data. This is useful when comparing results across different\n",
    "light curves on the same grid, or when you want fine control over the range.\n",
    "\n",
    "For `fast=True` the grid must be `np.linspace(0.0, max_omega, 2**k + 1)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": [
    "import light_curve as licu\n",
    "import numpy as np\n",
    "\n",
    "rng = np.random.default_rng(0)\n",
    "PERIOD = 23.5\n",
    "t = np.sort(rng.uniform(0, 300, 150))\n",
    "m = 15.0 + 0.4 * np.sin(2 * np.pi * t / PERIOD) + rng.normal(0, 0.04, 150)\n",
    "\n",
    "# Explicit grid: 0 to 2 rad/day in 2^10 + 1 steps (required for fast=True)\n",
    "omega_max = 2.0  # rad/day\n",
    "k = 10\n",
    "fixed_freqs = np.linspace(0.0, omega_max, 2 ** k + 1)\n",
    "\n",
    "pg_fixed = licu.Periodogram(peaks=1, freqs=fixed_freqs)\n",
    "result_fixed = pg_fixed(t, m)\n",
    "print(f'Best period (fixed grid): {result_fixed[0]:.3f} d')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16",
   "metadata": {},
   "source": "## Multiband periodogram\n\nWhen the same periodic signal is observed in several passbands, a joint\nLomb–Scargle periodogram pools all bands into a single frequency grid.\nBands with low amplitude still contribute, lifting the combined S/N above any\nsingle band.\n\nPass `bands=[...]` to `Periodogram` to switch to multiband mode.\nThe optional `multiband_normalization` (`'count'` or `'chi2'`) controls\nhow per-band powers are combined.\n\nThe example below synthesises a three-band light curve with two embedded\nperiods (P₁ = 23.5 d and P₂ = 7.3 d). Each band has different amplitudes,\nmagnitude offset, noise level, and a small inter-band phase shift (< 0.1 of the period):\n\n| band | amp(P₁) | amp(P₂) | offset | noise | phase(P₁) | phase(P₂) |\n|------|---------|---------|--------|-------|-----------|-----------|\n| g    | 0.35    | 0.10    | 15.0   | 0.04  | 0.00      | 0.00      |\n| r    | 0.22    | 0.18    | 15.5   | 0.15  | 0.03      | 0.05      |\n| i    | 0.12    | 0.28    | 15.9   | 0.06  | 0.07      | 0.08      |\n\nThe g band is dominated by P₁; the i band shows a stronger P₂ signal.\nThe r band is noisy (noise > amp(P₂)), so it cannot recover P₂ on its own.\nNo single band reliably recovers both periods, but the multiband periodogram does."
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "import light_curve as licu\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "BANDS = [\"g\", \"r\", \"i\"]\n",
    "COLORS = {\"g\": \"tab:green\", \"r\": \"tab:red\", \"i\": \"tab:purple\"}\n",
    "rng = np.random.default_rng(7)\n",
    "P1, P2 = 23.5, 7.3\n",
    "\n",
    "amp1 = {\"g\": 0.35, \"r\": 0.22, \"i\": 0.12}\n",
    "amp2 = {\"g\": 0.10, \"r\": 0.18, \"i\": 0.28}\n",
    "offset = {\"g\": 14.9, \"r\": 15.5, \"i\": 16.3}\n",
    "phi1 = {\"g\": 0.00, \"r\": 0.03, \"i\": 0.07}\n",
    "phi2 = {\"g\": 0.00, \"r\": 0.05, \"i\": 0.08}\n",
    "noise = {\"g\": 0.04, \"r\": 0.15, \"i\": 0.06}\n",
    "\n",
    "n_per_band = 80\n",
    "t_list, m_list, s_list, b_list = [], [], [], []\n",
    "for b in BANDS:\n",
    "    t_b = np.sort(rng.uniform(0, 400, n_per_band))\n",
    "    m_b = (offset[b]\n",
    "           + amp1[b] * np.sin(2 * np.pi * (t_b / P1 - phi1[b]))\n",
    "           + amp2[b] * np.sin(2 * np.pi * (t_b / P2 - phi2[b]))\n",
    "           + rng.normal(0, noise[b], n_per_band))\n",
    "    t_list.append(t_b)\n",
    "    m_list.append(m_b)\n",
    "    s_list.append(np.full(n_per_band, noise[b]))\n",
    "    b_list.append(np.full(n_per_band, b))\n",
    "\n",
    "t = np.concatenate(t_list)\n",
    "m = np.concatenate(m_list)\n",
    "sigma = np.concatenate(s_list)\n",
    "band = np.concatenate(b_list)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "for b in BANDS:\n",
    "    mask = band == b\n",
    "    phase = (t[mask] % P1) / P1\n",
    "    ax.errorbar(phase, m[mask], yerr=sigma[mask],\n",
    "                fmt='.', color=COLORS[b], alpha=0.6, ms=4, label=b)\n",
    "ax.invert_yaxis()\n",
    "ax.set_xlim(0, 1)\n",
    "ax.set_xlabel('Phase')\n",
    "ax.set_ylabel('mag')\n",
    "ax.set_title(f'Phase-folded at P₁ = {P1} d')\n",
    "ax.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18",
   "metadata": {},
   "outputs": [],
   "source": "pg_sb = licu.Periodogram(peaks=2, resolution=10, max_freq_factor=2.0, fast=True)\npg_mb = licu.Periodogram(peaks=2, bands=BANDS, resolution=10, max_freq_factor=2.0, fast=True)\n\nfig, axes = plt.subplots(4, 1, figsize=(10, 7), sharex=True)\n\nfor ax, b in zip(axes[:3], BANDS):\n    mask = band == b\n    omega, power = pg_sb.freq_power(t[mask], m[mask])\n    pos = omega > 0\n    periods = 2 * np.pi / omega[pos]\n    ax.plot(periods, power[pos], lw=0.7, color=COLORS[b], label=b)\n    for p_true in [P1, P2]:\n        ax.axvline(p_true, color='k', lw=0.8, ls='--', alpha=0.4)\n    ax.legend(loc='upper right')\n    ax.set_xscale('log')\n    ax.set_xlim(3, 60)\n\nomega_mb, power_mb = pg_mb.freq_power(t, m, sigma, band)\npos = omega_mb > 0\nperiods_mb = 2 * np.pi / omega_mb[pos]\naxes[3].plot(periods_mb, power_mb[pos], lw=0.7, color='k', label='multiband')\nfor p_true in [P1, P2]:\n    axes[3].axvline(p_true, color='k', lw=0.8, ls='--', alpha=0.4)\naxes[3].legend(loc='upper right')\naxes[3].set_xlabel('Period (days)')\nfig.suptitle('Lomb–Scargle power spectra  (dashed lines: true periods)')\nplt.tight_layout()\nplt.show()\n\nprint(f\"{'':>12}  {'P₀ (d)':>8}  {'S/N₀':>6}  {'P₁ (d)':>8}  {'S/N₁':>6}\")\nprint(\"-\" * 50)\nfor b in BANDS:\n    mask = band == b\n    r = pg_sb(t[mask], m[mask])\n    print(f\"{b:>12}  {r[0]:>8.2f}  {r[1]:>6.2f}  {r[2]:>8.2f}  {r[3]:>6.2f}\")\nr_mb = pg_mb(t, m, sigma, band=band)\nprint(\"-\" * 50)\nprint(f\"{'multiband':>12}  {r_mb[0]:>8.2f}  {r_mb[1]:>6.2f}  {r_mb[2]:>8.2f}  {r_mb[3]:>6.2f}\")\nprint(f\"{'(true)':>12}  {P1:>8.1f}  {'—':>6}  {P2:>8.1f}  {'—':>6}\")"
  },
  {
   "cell_type": "markdown",
   "id": "19",
   "metadata": {},
   "source": [
    "## See also\n",
    "\n",
    "- [Feature table](../) — full feature list\n",
    "- [Periodogram API](../api/periodogram/) — all constructor parameters"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}