{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.13.15","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":93783,"databundleVersionId":11331308,"sourceType":"competition"}],"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false},"papermill":{"default_parameters":{},"duration":14.216835,"end_time":"2025-03-10T09:42:54.211228","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2025-03-10T09:42:39.994393","version":"2.6.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"71eed97a","cell_type":"markdown","source":"# ESA-ADB Challenge — DAMP-style anomaly detection with STUMPY\n\nThis notebook uses a **left-looking matrix profile novelty detector**, using `stumpy.mass` to compute the z-normalized nearest-neighbor distance to earlier subsequences. This is the **DAMP score**.\n\nTo keep long ESA-ADB channels tractable, the defaults sample queries with `QUERY_STRIDE` and limit reference history with `MAX_HISTORY`. Set `QUERY_STRIDE=1` and `MAX_HISTORY=None` for exhaustive, uncapped left-profile scores (potentially *very slow* on this challenge). ","metadata":{}},{"id":"5e55d1a8","cell_type":"code","source":"# Install in Kaggle if needed\n!pip install -q stumpy","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-10-09T22:26:00.863741Z","iopub.execute_input":"2026-10-09T22:26:00.864007Z","iopub.status.idle":"2026-10-09T22:26:06.418256Z","shell.execute_reply.started":"2026-10-09T22:26:00.863974Z","shell.execute_reply":"2026-10-09T22:26:06.416994Z"}},"outputs":[],"execution_count":null},{"id":"984c7441","cell_type":"code","source":"import os\nimport time\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport stumpy\n\nprint('STUMPY:', stumpy.__version__)\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-10-09T22:26:06.42084Z","iopub.execute_input":"2026-10-09T22:26:06.421793Z","iopub.status.idle":"2026-10-09T22:26:10.554137Z","shell.execute_reply.started":"2026-10-09T22:26:06.421755Z","shell.execute_reply":"2026-10-09T22:26:10.553068Z"}},"outputs":[],"execution_count":null},{"id":"50891f2f","cell_type":"markdown","source":"## 1. Input data and parameters\n\n- `WINDOW_SIZE`: subsequence length in samples\n- `QUERY_STRIDE`: score one query every N samples; 1 is exhaustive.\n- `MAX_HISTORY`: maximum number of preceding reference samples (`None` = all). Larger values are slower.\n- `EXCLUSION_ZONE`: at least `WINDOW_SIZE // 4` preceding starts are excluded to prevent trivial overlap matches; adjacent windows should not be matched as their own reference.\n- `ANOMALY_QUANTILE`: unsupervised upper-tail threshold among valid computed query scores.\n- `EXPAND_WINDOWS`: map anomalous subsequences to all their covered timestamps instead of only their final sample.","metadata":{}},{"id":"b582080f","cell_type":"code","source":"CHANNEL_NAMES = ['channel_41']\nWINDOW_SIZE = 128\nQUERY_STRIDE = 16  # 1 computes every valid window (much slower)\nMAX_HISTORY = 10_000  # None uses the full preceding test history (very slow)\nEXCLUSION_ZONE = max(1, WINDOW_SIZE // 4)\nANOMALY_QUANTILE = 0.995\nEXPAND_WINDOWS = True\n\ncolumns_to_load = ['id'] + CHANNEL_NAMES + ['is_anomaly']\ntrain_df = pd.read_parquet('/kaggle/input/esa-adb-challenge/train.parquet', columns=columns_to_load)\ntest_df = pd.read_parquet('/kaggle/input/esa-adb-challenge/test.parquet', columns=['id'] + CHANNEL_NAMES)\nprint('Train:', train_df.shape, 'Test:', test_df.shape)\ndisplay(test_df.head())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-10-09T22:26:10.555319Z","iopub.execute_input":"2026-10-09T22:26:10.555955Z","iopub.status.idle":"2026-10-09T22:26:13.142525Z","shell.execute_reply.started":"2026-10-09T22:26:10.555916Z","shell.execute_reply":"2026-10-09T22:26:13.141697Z"}},"outputs":[],"execution_count":null},{"id":"70a83539","cell_type":"markdown","source":"## 2. Compute DAMP-style left matrix profile\n\nFor every sampled query window starting at index `i`, find its closest nontrivial **earlier** subsequence using `stumpy.mass`. For query `T[i:i+m]`, allowed reference starts end no later than `i - EXCLUSION_ZONE`. The earliest windows with insufficient left history have no score.","metadata":{}},{"id":"4ae59510","cell_type":"code","source":"def clean_channel(values):\n    series = pd.Series(np.asarray(values, dtype=np.float64))\n    series = series.replace([np.inf, -np.inf], np.nan)\n    if not series.notna().any():\n        raise ValueError('Selected channel has no finite observations')\n    return series.interpolate(limit_direction='both').to_numpy(dtype=np.float64)\n\n\ndef left_damp_profile(T, m, stride=1, max_history=None, exclusion_zone=None):\n    \"\"\"Compute sampled DAMP-style left-profile using STUMPY MASS.\n\n    Output indices identify *window starts*; NaNs indicate unevaluated positions.\n    For stride=1 and unlimited history, this is the exhaustive left-profile\n    distance (without the pruning optimizations of the published DAMP algorithm).\n    \"\"\"\n    T = np.asarray(T, dtype=np.float64)\n    n = len(T)\n    if n < 2 * m:\n        raise ValueError(f'Need at least 2*m samples; n={n}, m={m}')\n    ez = max(1, m // 4) if exclusion_zone is None else int(exclusion_zone)\n    scores = np.full(n - m + 1, np.nan, dtype=np.float64)\n    first_i = m + ez - 1  # reference window starting at zero is now eligible\n    for i in range(first_i, n - m + 1, stride):\n        latest_ref_start = i - ez\n        ref_start = 0 if max_history is None else max(0, i - max_history)\n        # Reference T[ref_start:latest_ref_start+m] includes candidate starts\n        # ref_start ... latest_ref_start (inclusive).\n        reference = T[ref_start:latest_ref_start + m]\n        if len(reference) < m:\n            continue\n        distances = stumpy.mass(T[i:i + m], reference)\n        distances = np.asarray(distances, dtype=np.float64)\n        finite = distances[np.isfinite(distances)]\n        if finite.size:\n            scores[i] = finite.min()\n    return scores\n\nstart = time.perf_counter()\nprofiles = {}\nfor channel in CHANNEL_NAMES:\n    print(f'Computing left matrix profile: {channel}...', flush=True)\n    profiles[channel] = left_damp_profile(\n        clean_channel(test_df[channel].to_numpy()), WINDOW_SIZE,\n        stride=QUERY_STRIDE, max_history=MAX_HISTORY,\n        exclusion_zone=EXCLUSION_ZONE,\n    )\nprint(f'Finished in {(time.perf_counter()-start)/60:.1f} minutes')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-10-09T22:26:13.143673Z","iopub.execute_input":"2026-10-09T22:26:13.144051Z","iopub.status.idle":"2026-10-09T22:28:37.045403Z","shell.execute_reply.started":"2026-10-09T22:26:13.144025Z","shell.execute_reply":"2026-10-09T22:28:37.044674Z"}},"outputs":[],"execution_count":null},{"id":"ad7168cd","cell_type":"markdown","source":"## 3. Convert subsequence scores into pointwise binary detections\n\nA high DAMP score signals a subsequence with no sufficiently similar earlier pattern. Each channel is thresholded by its own valid-score quantile, and detections are combined with a logical OR. When expanding a detected window to its covered points, use a difference array for efficiency. Early samples with no left history remain unflagged. The quantile is computed on the test series; do not use this setup for prospective/online detection without changing the threshold calibration.","metadata":{}},{"id":"e465f236","cell_type":"code","source":"n_samples = len(test_df)\ncombined_detections = np.zeros(n_samples, dtype=bool)\nthresholds = {}\nwindow_detections = {}\n\nfor channel, profile in profiles.items():\n    valid = np.isfinite(profile)\n    if not valid.any():\n        raise RuntimeError(f'No valid DAMP scores for {channel}; reduce WINDOW_SIZE or EXCLUSION_ZONE')\n    threshold = np.quantile(profile[valid], ANOMALY_QUANTILE)\n    thresholds[channel] = float(threshold)\n    starts = np.flatnonzero(valid & (profile > threshold))\n    window_detections[channel] = starts\n    if EXPAND_WINDOWS:\n        difference = np.zeros(n_samples + 1, dtype=np.int64)\n        np.add.at(difference, starts, 1)\n        np.add.at(difference, starts + WINDOW_SIZE, -1)\n        channel_flags = np.cumsum(difference[:-1]) > 0\n    else:\n        channel_flags = np.zeros(n_samples, dtype=bool)\n        channel_flags[starts + WINDOW_SIZE - 1] = True\n    combined_detections |= channel_flags\n    print(f'{channel}: threshold={threshold:.4f}, '\n          f'scored windows={valid.sum():,}, detected windows={len(starts):,}')\n\nsubmission_df = pd.DataFrame({\n    'id': test_df['id'].to_numpy(),\n    'is_anomaly': combined_detections.astype(np.uint8),\n})\nprint('Flagged points:', int(submission_df['is_anomaly'].sum()), '/', n_samples)\ndisplay(submission_df.head())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-10-09T22:28:37.046506Z","iopub.execute_input":"2026-10-09T22:28:37.046819Z","iopub.status.idle":"2026-10-09T22:28:37.095394Z","shell.execute_reply.started":"2026-10-09T22:28:37.046793Z","shell.execute_reply":"2026-10-09T22:28:37.094674Z"}},"outputs":[],"execution_count":null},{"id":"035fd866","cell_type":"markdown","source":"## 4. Write Kaggle submission","metadata":{}},{"id":"350fd67e","cell_type":"code","source":"os.makedirs('/kaggle/working', exist_ok=True)\nsubmission_df.to_csv('/kaggle/working/submission.csv', index=False)\nsubmission_df.to_parquet('/kaggle/working/DAMP-channel41.parquet', index=False)\nprint('Saved /kaggle/working/submission.csv')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-10-09T22:28:37.097261Z","iopub.execute_input":"2026-10-09T22:28:37.097547Z","iopub.status.idle":"2026-10-09T22:28:37.657274Z","shell.execute_reply.started":"2026-10-09T22:28:37.097523Z","shell.execute_reply":"2026-10-09T22:28:37.656427Z"}},"outputs":[],"execution_count":null},{"id":"37bf9ddf","cell_type":"markdown","source":"## 5. Inspect anomaly scores and detections","metadata":{}},{"id":"f59984d0","cell_type":"code","source":"channel = CHANNEL_NAMES[0]\nprofile = profiles[channel]\nwindow_ends = np.arange(len(profile)) + WINDOW_SIZE - 1\nvalid = np.isfinite(profile)\n\nfig, axes = plt.subplots(3, 1, figsize=(15, 9), sharex=True)\naxes[0].plot(np.arange(n_samples), test_df[channel].to_numpy(), lw=0.6)\naxes[0].set_ylabel(channel)\naxes[0].set_title('ESA-ADB test channel and DAMP-style novelty scores')\naxes[1].plot(window_ends[valid], profile[valid], lw=0.75, label='DAMP-style left profile')\naxes[1].axhline(thresholds[channel], ls='--', color='tab:red', label='Quantile threshold')\naxes[1].set_ylabel('Distance')\naxes[1].legend()\naxes[2].plot(np.arange(n_samples), submission_df['is_anomaly'].to_numpy(), lw=0.8)\naxes[2].set_ylabel('is_anomaly')\naxes[2].set_xlabel('Test sample index')\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-10-09T22:28:37.658292Z","iopub.execute_input":"2026-10-09T22:28:37.658534Z","iopub.status.idle":"2026-10-09T22:28:38.497178Z","shell.execute_reply.started":"2026-10-09T22:28:37.65851Z","shell.execute_reply":"2026-10-09T22:28:38.496294Z"}},"outputs":[],"execution_count":null},{"id":"33ef536a-e796-41cc-b3da-53a2e42ffee9","cell_type":"code","source":"fig, axes = plt.subplots(3, 1, figsize=(15, 9), sharex=True)\naxes[0].plot(np.arange(n_samples), test_df[channel].to_numpy(), lw=0.6)\naxes[0].set_ylabel(channel)\naxes[0].set_title('ESA-ADB test channel and DAMP-style novelty scores')\naxes[1].plot(window_ends[valid], profile[valid], lw=0.75, label='DAMP-style left profile')\naxes[1].axhline(thresholds[channel], ls='--', color='tab:red', label='Quantile threshold')\naxes[1].set_ylabel('Distance')\naxes[1].legend()\naxes[2].plot(np.arange(n_samples), submission_df['is_anomaly'].to_numpy(), lw=0.8)\naxes[2].set_ylabel('is_anomaly')\naxes[2].set_xlabel('Test sample index')\nplt.xlim([150000, 170000])\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-10-09T22:45:17.797496Z","iopub.execute_input":"2026-10-09T22:45:17.797865Z","iopub.status.idle":"2026-10-09T22:45:18.341737Z","shell.execute_reply.started":"2026-10-09T22:45:17.797838Z","shell.execute_reply":"2026-10-09T22:45:18.340744Z"}},"outputs":[],"execution_count":null}]}