{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "56a0cc6c",
   "metadata": {},
   "source": [
    "# MaldiAMRKit - Peak Sets\n",
    "\n",
    "This notebook covers the **peak sets**: compact `(m/z, intensity)` containers (`PeakSet` / `PeakList`) that represent each spectrum as a variable-length set of peaks for downstream peak-set models. It shows how to extract peak sets, rank and pad peaks, persist and cache them, and apply fit-free peak alignment."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b4f51f2",
   "metadata": {},
   "source": [
    "## Import Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7a5d2bb9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-25T09:47:51.191304Z",
     "iopub.status.busy": "2026-06-25T09:47:51.191209Z",
     "iopub.status.idle": "2026-06-25T09:47:52.011174Z",
     "shell.execute_reply": "2026-06-25T09:47:52.010598Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from maldiamrkit import MaldiSet\n",
    "from maldiamrkit.alignment import align_peaks\n",
    "from maldiamrkit.detection import MaldiPeakDetector, PeakList, create_peakset_input"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bd600e9",
   "metadata": {},
   "source": [
    "## Load Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "66ace42b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-25T09:47:52.012869Z",
     "iopub.status.busy": "2026-06-25T09:47:52.012701Z",
     "iopub.status.idle": "2026-06-25T09:47:52.622129Z",
     "shell.execute_reply": "2026-06-25T09:47:52.621605Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Binned spectra shape: (29, 6000)\n"
     ]
    }
   ],
   "source": [
    "data = MaldiSet.from_directory(\n",
    "    \"../data/\",\n",
    "    \"../data/metadata/metadata.csv\",\n",
    "    aggregate_by=dict(antibiotics=\"Drug\"),\n",
    ")\n",
    "X = data.X\n",
    "print(f\"Binned spectra shape: {X.shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99fc33f4",
   "metadata": {},
   "source": [
    "## Extract a PeakList from binned spectra\n",
    "\n",
    "`MaldiPeakDetector.transform_peaklist` detects peaks per spectrum and returns a `PeakList` - one `PeakSet` per spectrum, each a set of `(m/z, intensity)` peaks. Extraction is **stateless** (a pure per-spectrum function), so the result is identical whether computed over the whole dataset or inside a single cross-validation fold. The m/z axis is recovered from the binned column labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ec0c6710",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-25T09:47:52.623274Z",
     "iopub.status.busy": "2026-06-25T09:47:52.623175Z",
     "iopub.status.idle": "2026-06-25T09:47:52.633910Z",
     "shell.execute_reply": "2026-06-25T09:47:52.633465Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PeakList(n=29, method=local, top_k=50)\n",
      "Spectra: 29\n",
      "Per-spectrum peak counts: [50 50 50 50 50 50 50 50] ...\n",
      "warped (pre-alignment): False\n"
     ]
    }
   ],
   "source": [
    "detector = MaldiPeakDetector(method=\"local\", prominence=1e-4)\n",
    "peaks = detector.transform_peaklist(X, top_k=50)\n",
    "\n",
    "print(peaks)\n",
    "print(f\"Spectra: {len(peaks)}\")\n",
    "print(f\"Per-spectrum peak counts: {peaks.n_peaks[:8]} ...\")\n",
    "print(f\"warped (pre-alignment): {peaks.meta['warped']}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b9a33c3",
   "metadata": {},
   "source": [
    "## Inspect a single PeakSet\n",
    "\n",
    "Peaks are stored sorted by ascending m/z. Index the `PeakList` by position to get a `PeakSet`; `as_array()` returns an `(n_peaks, 2)` array of `(m/z, intensity)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "71d2137d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-25T09:47:52.634848Z",
     "iopub.status.busy": "2026-06-25T09:47:52.634763Z",
     "iopub.status.idle": "2026-06-25T09:47:52.636944Z",
     "shell.execute_reply": "2026-06-25T09:47:52.636429Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PeakSet(n_peaks=50, mz=[2012.0..9530.0])\n",
      "n_peaks: 50\n",
      "first 5 (m/z, intensity):\n",
      "[[2.01200000e+03 6.87894941e-04]\n",
      " [2.03000000e+03 3.54902016e-03]\n",
      " [2.04500000e+03 9.10884239e-04]\n",
      " [2.05100000e+03 1.97589893e-03]\n",
      " [2.06000000e+03 9.63189503e-04]]\n"
     ]
    }
   ],
   "source": [
    "ps = peaks[0]\n",
    "print(ps)\n",
    "print(f\"n_peaks: {ps.n_peaks}\")\n",
    "print(f\"first 5 (m/z, intensity):\\n{ps.as_array()[:5]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd5e813d",
   "metadata": {},
   "source": [
    "## Rank and keep the strongest peaks\n",
    "\n",
    "`top_k` keeps the highest-ranked peaks. The ranking is selected at detection time via `rank_by`: `\"intensity\"` (default), `\"prominence\"`, or `\"persistence\"`. For a non-intensity ranking, each `PeakSet` carries the per-peak `score`, so any later truncation (e.g. padding) stays consistent with how peaks were detected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a7e9ed51",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-25T09:47:52.637836Z",
     "iopub.status.busy": "2026-06-25T09:47:52.637738Z",
     "iopub.status.idle": "2026-06-25T09:47:52.648724Z",
     "shell.execute_reply": "2026-06-25T09:47:52.648189Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rank_by: prominence\n",
      "per-peak score carried: True\n",
      "5 strongest m/z in spectrum 0: [2030. 2201. 4361. 5375. 6284.]\n"
     ]
    }
   ],
   "source": [
    "peaks_prom = detector.transform_peaklist(X, top_k=20, rank_by=\"prominence\")\n",
    "ps_prom = peaks_prom[0]\n",
    "\n",
    "print(f\"rank_by: {peaks_prom.meta['rank_by']}\")\n",
    "print(f\"per-peak score carried: {ps_prom.score is not None}\")\n",
    "print(f\"5 strongest m/z in spectrum 0: {ps_prom.top_k(5).mz}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d09e718",
   "metadata": {},
   "source": [
    "## Persistent-homology ranking\n",
    "\n",
    "Persistent homology (`method=\"ph\"`) ranks peaks by topological persistence, which is robust to noise. Persistence ranking requires the PH detector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a9acbf09",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-25T09:47:52.649682Z",
     "iopub.status.busy": "2026-06-25T09:47:52.649597Z",
     "iopub.status.idle": "2026-06-25T09:47:52.822549Z",
     "shell.execute_reply": "2026-06-25T09:47:52.821931Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PeakList(n=29, method=ph, top_k=20)\n",
      "persistence scores (spectrum 0, first 5): [0.00354396 0.00167008 0.00120384 0.00132015 0.0027002 ]\n"
     ]
    }
   ],
   "source": [
    "ph_detector = MaldiPeakDetector(method=\"ph\", persistence_threshold=5e-4)\n",
    "peaks_ph = ph_detector.transform_peaklist(X, top_k=20, rank_by=\"persistence\")\n",
    "\n",
    "print(peaks_ph)\n",
    "print(f\"persistence scores (spectrum 0, first 5): {peaks_ph[0].score[:5]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c254a6a",
   "metadata": {},
   "source": [
    "## Pad to a dense array for modelling\n",
    "\n",
    "Peak-set models consume a fixed-width tensor. `to_padded` packs the variable-length sets into a dense `(n_spectra, P, 2)` array plus a boolean mask flagging real (non-padded) peaks. Spectra with more than `P` peaks keep their top-ranked ones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "28205869",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-25T09:47:52.823783Z",
     "iopub.status.busy": "2026-06-25T09:47:52.823693Z",
     "iopub.status.idle": "2026-06-25T09:47:52.825977Z",
     "shell.execute_reply": "2026-06-25T09:47:52.825463Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "values shape: (29, 50, 2)  (n_spectra, P, 2 = m/z, intensity)\n",
      "mask shape:   (29, 50)\n",
      "real peaks in spectrum 0: 50 / 50\n"
     ]
    }
   ],
   "source": [
    "values, mask = peaks.to_padded(max_peaks=50)\n",
    "\n",
    "print(f\"values shape: {values.shape}  (n_spectra, P, 2 = m/z, intensity)\")\n",
    "print(f\"mask shape:   {mask.shape}\")\n",
    "print(f\"real peaks in spectrum 0: {int(mask[0].sum())} / {mask.shape[1]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b06c5aff",
   "metadata": {},
   "source": [
    "## Index and slice a PeakList\n",
    "\n",
    "Integer access returns a `PeakSet`; slicing returns a sub-`PeakList` that keeps the matching sample index and metadata."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6f52855a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-25T09:47:52.826991Z",
     "iopub.status.busy": "2026-06-25T09:47:52.826910Z",
     "iopub.status.idle": "2026-06-25T09:47:52.828915Z",
     "shell.execute_reply": "2026-06-25T09:47:52.828418Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PeakList(n=5, method=local, top_k=50)\n",
      "index: ['10s', '11s', '12s', '13s', '14s']\n"
     ]
    }
   ],
   "source": [
    "subset = peaks[:5]\n",
    "print(subset)\n",
    "print(f\"index: {list(subset.index)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c637f7dd",
   "metadata": {},
   "source": [
    "## Persist and cache\n",
    "\n",
    "`save` / `load` round-trips a `PeakList` to `.npz` + `.json`. Because extraction is stateless and leak-safe, a `PeakList` can also be cached: pass `cache_dir=` and a repeat call with the same data and configuration reuses the stored result instead of recomputing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "559dba9b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-25T09:47:52.829774Z",
     "iopub.status.busy": "2026-06-25T09:47:52.829693Z",
     "iopub.status.idle": "2026-06-25T09:47:52.842211Z",
     "shell.execute_reply": "2026-06-25T09:47:52.841740Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reloaded 29 spectra (method=local)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cached files: ['peaklist_9593a4e7dea1d102_bdbc80b1d9a8bd5b.json', 'peaklist_9593a4e7dea1d102_bdbc80b1d9a8bd5b.npz']\n"
     ]
    }
   ],
   "source": [
    "import tempfile\n",
    "from pathlib import Path\n",
    "\n",
    "tmp = Path(tempfile.mkdtemp())\n",
    "\n",
    "# save / load round-trip\n",
    "peaks.save(tmp / \"peaks\")\n",
    "reloaded = PeakList.load(tmp / \"peaks\")\n",
    "print(f\"reloaded {len(reloaded)} spectra (method={reloaded.meta['method']})\")\n",
    "\n",
    "# content + config keyed cache: the second call is a cache hit\n",
    "cache = tmp / \"cache\"\n",
    "detector.transform_peaklist(X, top_k=50, cache_dir=cache)\n",
    "detector.transform_peaklist(X, top_k=50, cache_dir=cache)\n",
    "print(f\"cached files: {sorted(p.name for p in cache.glob('peaklist_*'))}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "712adb6d",
   "metadata": {},
   "source": [
    "## Binning-free extraction from raw files\n",
    "\n",
    "`create_peakset_input` is the faithful, binning-free path: it preprocesses each raw spectrum and detects peaks at full m/z resolution (no binning), returning an un-warped `PeakList`. It mirrors `create_raw_input` from the alignment module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "eb0b4d77",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-25T09:47:52.843562Z",
     "iopub.status.busy": "2026-06-25T09:47:52.843477Z",
     "iopub.status.idle": "2026-06-25T09:47:53.139497Z",
     "shell.execute_reply": "2026-06-25T09:47:53.138657Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PeakList(n=29, method=local, top_k=50)\n",
      "source: create_peakset_input\n",
      "warped: False\n"
     ]
    }
   ],
   "source": [
    "peaks_raw = create_peakset_input(\"../data/\", top_k=50)\n",
    "\n",
    "print(peaks_raw)\n",
    "print(f\"source: {peaks_raw.meta['source']}\")\n",
    "print(f\"warped: {peaks_raw.meta['warped']}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80ce4ef3",
   "metadata": {},
   "source": [
    "## Fit-free peak alignment\n",
    "\n",
    "`align_peaks` warps a peak set's m/z onto a reference peak list using the shared alignment strategies. It is **fit-free**: the caller supplies the reference (estimated on the training split, in-fold), so it never leaks. Intensities are carried through unchanged. Here we shift a spectrum onto a reference deliberately offset by 5 Da."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb1cf9d3",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-25T09:47:53.141106Z",
     "iopub.status.busy": "2026-06-25T09:47:53.141006Z",
     "iopub.status.idle": "2026-06-25T09:47:53.143393Z",
     "shell.execute_reply": "2026-06-25T09:47:53.143105Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "median m/z shift applied: 5.00 Da\n",
      "n_peaks unchanged: True\n",
      "intensities unchanged: True\n"
     ]
    }
   ],
   "source": [
    "sample = peaks_raw[0]\n",
    "ref_mz = sample.mz + 5.0  # generate an offset of 5 Da\n",
    "\n",
    "aligned = align_peaks(sample, ref_mz, method=\"shift\", max_shift_da=20.0)\n",
    "\n",
    "print(f\"median m/z shift applied: {np.median(aligned.mz - sample.mz):.2f} Da\")\n",
    "print(f\"n_peaks unchanged: {aligned.n_peaks == sample.n_peaks}\")\n",
    "print(f\"intensities unchanged: {np.array_equal(sample.intensity, aligned.intensity)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7fc0c41",
   "metadata": {},
   "source": [
    "## Avoiding data leakage\n",
    "\n",
    "Peak extraction is a pure per-spectrum function (`meta[\"warped\"] = False`), so precomputing or caching a `PeakList` over the whole dataset cannot leak. Any *fitted* step (an `align_peaks` reference, standardization) must still be estimated on the training split, inside the fold."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "maldiamrkit",
   "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}