{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# dm–dt Tutorial\n",
    "\n",
    "A dm-dt map is a 2D histogram of magnitude differences (Δm) versus log-time differences\n",
    "(lg Δt) for all pairs of observations in a light curve.\n",
    "It was introduced as an ML input representation by\n",
    "[Mahabal et al. 2011](https://ui.adsabs.harvard.edu/abs/2011BASI...39..387M/abstract).\n",
    "\n",
    "`DmDt` produces a fixed-size 2D array that can be fed directly into a CNN."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %pip install light-curve matplotlib"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import light_curve as licu\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as mcolors\n",
    "\n",
    "dmdt = licu.DmDt.from_borders(\n",
    "    min_lgdt=0,      # lg Δt_min = 1 day\n",
    "    max_lgdt=2,      # lg Δt_max = 100 days\n",
    "    max_abs_dm=1.5,  # |Δm| up to 1.5 mag\n",
    "    lgdt_size=32,\n",
    "    dm_size=32,\n",
    "    norm=[\"dt\"],\n",
    ")\n",
    "print(f\"Map shape: {dmdt.shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "## Computing a map\n",
    "\n",
    "Call `.points(t, m)` to compute the dm-dt map for a single light curve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(0)\n",
    "t = np.sort(rng.uniform(0, 100, 200)).astype(np.float64)\n",
    "\n",
    "period = 7.0\n",
    "m = 15.0 + 0.2 * np.sin(2 * np.pi * t / period) + rng.normal(0, 0.03, 200)\n",
    "\n",
    "map_ = dmdt.points(t, m)\n",
    "print(f\"Map shape: {map_.shape}  (lgdt_size × dm_size)\")\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5, 4))\n",
    "im = ax.imshow(\n",
    "    map_.T,\n",
    "    origin=\"lower\",\n",
    "    aspect=\"auto\",\n",
    "    extent=[0, 2, -1.5, 1.5],\n",
    "    cmap=\"Blues\",\n",
    "    norm=mcolors.PowerNorm(gamma=0.5, vmin=0),\n",
    ")\n",
    "ax.set_xlabel(\"lg Δt (days)\")\n",
    "ax.set_ylabel(\"Δm (mag)\")\n",
    "ax.set_title(\"dm-dt map — periodic variable (P = 7 d)\")\n",
    "fig.colorbar(im, ax=ax, label=\"normalised count\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "## Comparing light curve types\n",
    "\n",
    "Different variability classes leave distinct signatures in dm-dt space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng2 = np.random.default_rng(42)\n",
    "t2 = np.sort(rng2.uniform(0, 100, 200)).astype(np.float64)\n",
    "\n",
    "configs = [\n",
    "    (\"Constant star\",        15.0 + rng2.normal(0, 0.05, 200)),\n",
    "    (\"Periodic (P = 7 d)\",  15.0 + 0.4 * np.sin(2 * np.pi * t2 / 7) + rng2.normal(0, 0.03, 200)),\n",
    "    (\"Transient (Gaussian)\", 15.0 - 1.0 * np.exp(-((t2 - 50) ** 2) / 50) + rng2.normal(0, 0.05, 200)),\n",
    "]\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(13, 4), sharey=True)\n",
    "norm = mcolors.PowerNorm(gamma=0.5, vmin=0)\n",
    "\n",
    "for ax, (title, m2) in zip(axes, configs):\n",
    "    mp = dmdt.points(t2, m2)\n",
    "    im = ax.imshow(\n",
    "        mp.T, origin=\"lower\", aspect=\"auto\",\n",
    "        extent=[0, 2, -1.5, 1.5], cmap=\"Blues\", norm=norm,\n",
    "    )\n",
    "    ax.set_title(title)\n",
    "    ax.set_xlabel(\"lg Δt (days)\")\n",
    "    fig.colorbar(im, ax=ax, label=\"normalised count\", shrink=0.8)\n",
    "\n",
    "axes[0].set_ylabel(\"Δm (mag)\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "## Error-weighted maps (`gausses`)\n",
    "\n",
    "When photometric errors are known, use `.gausses()` to spread each observation pair\n",
    "into a Gaussian kernel in Δm space. The third argument is the **variance** σ² (not σ).\n",
    "\n",
    "Here we use strong noise (σ = 0.3 mag) to make the difference clearly visible:\n",
    "\n",
    "!!! warning \"Variance, not std dev\"\n",
    "     expects the **variance** σ² — remember to square your error array: .\n",
    "\n",
    "!!! warning \"Variance, not std dev\"\n",
    "    `.gausses()` expects the **variance** σ² — remember to square your error array: `err ** 2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng3 = np.random.default_rng(7)\n",
    "t3 = np.sort(rng3.uniform(0, 100, 150)).astype(np.float64)\n",
    "err_large = 0.3  # large photometric error to make the effect visible\n",
    "m3 = 15.0 + 0.5 * np.sin(2 * np.pi * t3 / 7) + rng3.normal(0, err_large, 150)\n",
    "err3 = np.full(150, err_large)\n",
    "\n",
    "map_pts  = dmdt.points(t3, m3)\n",
    "map_gaus = dmdt.gausses(t3, m3, err3 ** 2)  # err² = variance\n",
    "map_diff = map_pts - map_gaus\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(13, 4), sharey=True)\n",
    "kw = dict(origin=\"lower\", aspect=\"auto\", extent=[0, 2, -1.5, 1.5],\n",
    "          cmap=\"Blues\", norm=mcolors.PowerNorm(gamma=0.5, vmin=0))\n",
    "\n",
    "im1 = axes[0].imshow(map_pts.T,  **kw)\n",
    "axes[0].set_title(\"points()\")\n",
    "fig.colorbar(im1, ax=axes[0], label=\"count\", shrink=0.8)\n",
    "\n",
    "im2 = axes[1].imshow(map_gaus.T, **kw)\n",
    "axes[1].set_title(\"gausses() — error-weighted\")\n",
    "fig.colorbar(im2, ax=axes[1], label=\"count\", shrink=0.8)\n",
    "\n",
    "# Difference panel: use a diverging colormap\n",
    "vmax = np.abs(map_diff).max()\n",
    "im3 = axes[2].imshow(map_diff.T, origin=\"lower\", aspect=\"auto\",\n",
    "                     extent=[0, 2, -1.5, 1.5],\n",
    "                     cmap=\"RdBu_r\", vmin=-vmax, vmax=vmax)\n",
    "axes[2].set_title(\"points() − gausses()\")\n",
    "fig.colorbar(im3, ax=axes[2], label=\"difference\", shrink=0.8)\n",
    "\n",
    "for ax in axes:\n",
    "    ax.set_xlabel(\"lg Δt (days)\")\n",
    "axes[0].set_ylabel(\"Δm (mag)\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "## Batch processing\n",
    "\n",
    "`.points_many()` computes maps for a list of light curves in one call,\n",
    "returning a 3D array of shape `(N, lgdt_size, dm_size)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng4 = np.random.default_rng(99)\n",
    "light_curves = []\n",
    "for _ in range(100):\n",
    "    n_obs = int(rng4.integers(50, 200))\n",
    "    t_i = np.sort(rng4.uniform(0, 100, n_obs)).astype(np.float64)\n",
    "    m_i = rng4.normal(15, 0.3, n_obs)\n",
    "    light_curves.append((t_i, m_i))\n",
    "\n",
    "maps = dmdt.points_many(light_curves)\n",
    "print(f\"Batch output shape: {maps.shape}  # (N, lgdt_size, dm_size)\")\n",
    "print(f\"Ready for CNNs or flatten to (N, {maps.shape[1] * maps.shape[2]}) for sklearn.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "metadata": {},
   "source": [
    "## Normalisation\n",
    "\n",
    "`norm` controls how the raw pair counts are scaled before returning the map.\n",
    "Three modes are available (combinable):\n",
    "\n",
    "| `norm` value | Effect |\n",
    "|---|---|\n",
    "| `[]` (default) | Raw pair counts |\n",
    "| `[\"dt\"]` | Each lg-Δt row divided by its total count — gives a **conditional distribution** p(Δm \\| Δt) per row |\n",
    "| `[\"max\"]` | Entire map divided by its maximum value — brings scale to [0, 1] |\n",
    "| `[\"dt\", \"max\"]` | `\"dt\"` first, then `\"max\"` applied to the result |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import light_curve as licu\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as mcolors\n",
    "\n",
    "rng = np.random.default_rng(0)\n",
    "t = np.sort(rng.uniform(0, 100, 200)).astype(np.float64)\n",
    "m = 15.0 + 0.5 * np.sin(2 * np.pi * t / 7.0) + rng.normal(0, 0.03, 200)\n",
    "\n",
    "norm_opts = [[], [\"dt\"], [\"max\"], [\"dt\", \"max\"]]\n",
    "titles = [\"no norm\", \"norm=[dt]\", \"norm=[max]\", \"norm=[dt,max]\"]\n",
    "\n",
    "fig, axes = plt.subplots(1, 4, figsize=(14, 4), sharey=True)\n",
    "\n",
    "for ax, norm, title in zip(axes, norm_opts, titles):\n",
    "    d = licu.DmDt.from_borders(\n",
    "        min_lgdt=0, max_lgdt=2, max_abs_dm=1.5,\n",
    "        lgdt_size=32, dm_size=32, norm=norm,\n",
    "    )\n",
    "    mp = d.points(t, m)\n",
    "    per_norm = mcolors.PowerNorm(gamma=0.5, vmin=0, vmax=mp.max() or 1)\n",
    "    im = ax.imshow(\n",
    "        mp.T, origin=\"lower\", aspect=\"auto\",\n",
    "        extent=[0, 2, -1.5, 1.5], cmap=\"Blues\", norm=per_norm,\n",
    "    )\n",
    "    fig.colorbar(im, ax=ax, shrink=0.8)\n",
    "    ax.set_xlabel(\"lg Δt (days)\")\n",
    "    ax.set_title(title)\n",
    "\n",
    "axes[0].set_ylabel(\"Δm (mag)\")\n",
    "plt.suptitle(\"Effect of dm–dt map normalisation\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14",
   "metadata": {},
   "source": [
    "## See also\n",
    "\n",
    "- [API reference](../api/)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}