{"nbformat": 4, "nbformat_minor": 5, "metadata": {"kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"}, "language_info": {"name": "python", "version": "3.10.12"}}, "cells": [{"cell_type": "markdown", "metadata": {}, "source": "# Corporaci\u00f3n Favorita \u2014 Grocery Sales Forecasting\n\n**Retrospective LightGBM benchmark** \u00b7 Kaggle 2017\u20132018 competition (closed)\n\n## What this notebook produces\n\nA 16-day sales forecast (2017-08-16 \u2192 2017-08-31) for all 210,654 store\u2013item pairs in the competition test set, using a 72-feature LightGBM pipeline blended across three random seeds.\n\n## Fork-and-run: what you reproduce\n\n| | |\n|---|---|\n| **Private LB NWRMSLE** | 0.52840 |\n| **Runtime** (Kaggle 2\u00d7 vCPU CPU kernel, no GPU) | ~6 h |\n| **Peak RAM** | ~17.5 GB (30 GB budget; stays single-member sequential) |\n| **Submission rows** | 3,370,464 |\n\nNo GPU required. No internet required. The competition dataset is auto-attached via `competition_sources` in the kernel metadata.\n\n**Smoke test:** set `SMOKE_MODE = True` in Cell 2 for a ~3-min end-to-end run (n_estimators=5, 2013 cutoff, dev stores only) that exercises every cell and asserts the row count.\n\n## Honesty\n\n- **Retrospective benchmark.** Hyperparameters were chosen with foreknowledge of backtest direction; this is not a blind prospective study.\n- **Published solutions existed** at development time (top-5 writeups and kernels were publicly available).\n- **Promotion schedule is competition-provided.** The `onpromo` column in `test.csv` gives promotion status during the forecast horizon. Loading it with `allow_future=True` is permitted by competition design (see the submission cell). In a real forecasting setting this information would not be available.\n- **`is_earthquake` dead by construction.** The 2016-04-16 earthquake feature is identically zero for all 2017 training cutoffs; it was removed rather than fixed."}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "import os, glob\nfrom pathlib import Path\n\n# Discover input dir dynamically \u2014 Kaggle nests competition data differently\n# depending on whether the dataset is attached as a competition or user dataset.\n# The recursive glob finds train.csv.7z at any nesting depth.\nhits = glob.glob('/kaggle/input/**/train.csv.7z', recursive=True)\nassert len(hits) == 1, (\n    f\"Expected exactly one train.csv.7z, found {len(hits)}: {hits}. \"\n    \"Check that the Favorita competition dataset is attached to this notebook.\"\n)\nINPUT_DIR = Path(hits[0]).parent\nprint(f\"Discovered INPUT_DIR: {INPUT_DIR}\")\n\nEXPECTED_ARCHIVES = [\n    \"train.csv.7z\",\n    \"test.csv.7z\",\n    \"stores.csv.7z\",\n    \"items.csv.7z\",\n    \"holidays_events.csv.7z\",\n    \"oil.csv.7z\",\n    \"transactions.csv.7z\",\n    \"sample_submission.csv.7z\",\n]\n\nprint(\"\\nInput directory listing:\")\nfor f in sorted(INPUT_DIR.iterdir()):\n    print(f\"  {f.name:<35s}  {f.stat().st_size / 1e6:>8.1f} MB\")\n\nmissing = [f for f in EXPECTED_ARCHIVES if not (INPUT_DIR / f).exists()]\nassert not missing, (\n    f\"Missing {len(missing)} expected archive(s): {missing}\\n\"\n    \"Check that the Favorita competition dataset is attached to this notebook.\"\n)\nprint(f\"\\nAll {len(EXPECTED_ARCHIVES)} expected .7z archives confirmed present.\")"}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "import random, os, sys, time\nimport subprocess\nimport numpy as np\nimport pandas as pd\nimport psutil\nimport pyarrow as pa\nimport pyarrow.parquet as pq\n\nSMOKE_MODE = False   # Set True for a fast end-to-end smoke test (~2 min)\n\nif SMOKE_MODE:\n    print(\"=\" * 60)\n    print(\"  SMOKE MODE \u2014 2-chunk train cap, n_estimators=5\")\n    print(\"  Every cell runs; row-count assert is live.\")\n    print(\"=\" * 60)\n\nSEED = 42\nrandom.seed(SEED)\nnp.random.seed(SEED)\nt_notebook_start = time.time()\n\n# Point dataguard at the working data directory\nos.environ[\"FAVORITA_DATA_ROOT\"] = \"/kaggle/working/data\"\n\nos.makedirs(\"/kaggle/working/data\", exist_ok=True)\nos.makedirs(\"/kaggle/working/dataguard\", exist_ok=True)\nos.makedirs(\"/kaggle/working/models\", exist_ok=True)\n\nprint(\"Directories ready. SEED =\", SEED)"}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "from pathlib import Path\nimport time as _time\n\nRAW  = INPUT_DIR   # discovered in cell 0 via glob\nDATA = Path(\"/kaggle/working/data\")\nCHUNKSIZE = 500_000\n\nTRAIN_DTYPES = {\"id\": \"int32\", \"store_nbr\": \"int8\", \"item_nbr\": \"int32\", \"unit_sales\": \"float32\"}\nTEST_DTYPES  = {\"id\": \"int32\", \"store_nbr\": \"int8\", \"item_nbr\": \"int32\"}\nTRANSACTIONS_DTYPES = {\"store_nbr\": \"int8\", \"transactions\": \"int32\"}\nITEMS_DTYPES = {\"item_nbr\": \"int32\", \"class\": \"int16\", \"perishable\": \"int8\"}\n\n\ndef _7z_stdout(archive_name: str) -> subprocess.Popen:\n    \"\"\"Stream-decompress a .7z archive from the Kaggle input dir to stdout.\n    7z (p7zip-full) is preinstalled in the Kaggle Docker image \u2014 no Internet needed.\"\"\"\n    return subprocess.Popen(\n        [\"7z\", \"e\", str(RAW / archive_name), \"-so\"],\n        stdout=subprocess.PIPE,\n        stderr=subprocess.DEVNULL,\n    )\n\n\nt0 = _time.time()\n\n# \u2500\u2500 stores \u2500\u2500\nproc = _7z_stdout(\"stores.csv.7z\")\ndf = pd.read_csv(proc.stdout)\nproc.wait()\ndf[\"store_nbr\"] = df[\"store_nbr\"].astype(\"int8\")\ndf[\"cluster\"]   = df[\"cluster\"].astype(\"int8\")\ndf[\"type\"]      = df[\"type\"].astype(\"category\")\ndf[\"city\"]      = df[\"city\"].astype(\"category\")\ndf[\"state\"]     = df[\"state\"].astype(\"category\")\ndf.to_parquet(DATA / \"stores.parquet\", index=False)\nprint(f\"stores:       {len(df):,} rows\")\n\n# \u2500\u2500 items \u2500\u2500\nproc = _7z_stdout(\"items.csv.7z\")\ndf = pd.read_csv(proc.stdout, dtype=ITEMS_DTYPES)\nproc.wait()\ndf[\"family\"] = df[\"family\"].astype(\"category\")\ndf.to_parquet(DATA / \"items.parquet\", index=False)\nprint(f\"items:        {len(df):,} rows\")\n\n# \u2500\u2500 holidays \u2500\u2500\nproc = _7z_stdout(\"holidays_events.csv.7z\")\ndf = pd.read_csv(proc.stdout, parse_dates=[\"date\"])\nproc.wait()\ndf[\"type\"]        = df[\"type\"].astype(\"category\")\ndf[\"locale\"]      = df[\"locale\"].astype(\"category\")\ndf[\"locale_name\"] = df[\"locale_name\"].astype(\"category\")\ndf[\"transferred\"] = df[\"transferred\"].map({\"True\": True, \"False\": False}).astype(\"boolean\")\ndf.to_parquet(DATA / \"holidays.parquet\", index=False)\nprint(f\"holidays:     {len(df):,} rows\")\n\n# \u2500\u2500 oil \u2500\u2500\nproc = _7z_stdout(\"oil.csv.7z\")\ndf = pd.read_csv(proc.stdout, parse_dates=[\"date\"])\nproc.wait()\ndf[\"dcoilwtico\"] = df[\"dcoilwtico\"].astype(\"float32\")\ndf.to_parquet(DATA / \"oil.parquet\", index=False)\nprint(f\"oil:          {len(df):,} rows\")\n\n# \u2500\u2500 transactions \u2500\u2500\nproc = _7z_stdout(\"transactions.csv.7z\")\ndf = pd.read_csv(proc.stdout, parse_dates=[\"date\"], dtype=TRANSACTIONS_DTYPES)\nproc.wait()\ndf.to_parquet(DATA / \"transactions.parquet\", index=False)\nprint(f\"transactions: {len(df):,} rows\")\n\n# \u2500\u2500 test \u2500\u2500\nproc = _7z_stdout(\"test.csv.7z\")\ndf = pd.read_csv(proc.stdout, dtype=TEST_DTYPES, parse_dates=[\"date\"])\nproc.wait()\ndf[\"onpromotion\"] = df[\"onpromotion\"].astype(\"boolean\")\ndf.to_parquet(DATA / \"test.parquet\", index=False)\nprint(f\"test:         {len(df):,} rows\")\n\n# \u2500\u2500 sample_submission \u2500\u2500\nproc = _7z_stdout(\"sample_submission.csv.7z\")\ndf = pd.read_csv(proc.stdout, dtype={\"id\": \"int32\", \"unit_sales\": \"float32\"})\nproc.wait()\ndf.to_parquet(DATA / \"sample_submission.parquet\", index=False)\nprint(f\"sample_sub:   {len(df):,} rows\")\n\n# \u2500\u2500 train (chunked \u2014 ~125M rows full / 2 chunks in SMOKE_MODE) \u2500\u2500\nproc = _7z_stdout(\"train.csv.7z\")\nout_train  = DATA / \"train.parquet\"\nwriter     = None\ntotal_rows = 0\nt_train    = _time.time()\nreader = pd.read_csv(\n    proc.stdout,\n    dtype=TRAIN_DTYPES,\n    parse_dates=[\"date\"],\n    chunksize=CHUNKSIZE,\n    na_values=[\"\"],\n    keep_default_na=True,\n)\nfor i, chunk in enumerate(reader):\n    chunk[\"onpromotion\"] = chunk[\"onpromotion\"].astype(\"boolean\")\n    table = pa.Table.from_pandas(chunk, preserve_index=False)\n    if writer is None:\n        writer = pq.ParquetWriter(out_train, table.schema, compression=\"snappy\")\n    writer.write_table(table)\n    total_rows += len(chunk)\n    if (i + 1) % 50 == 0:\n        print(f\"  train chunk {i+1:4d}: {total_rows:,} rows  {_time.time()-t_train:.0f}s\", flush=True)\n    if SMOKE_MODE and (i + 1) >= 2:   # cap at 2 chunks (~1M rows) for smoke test\n        break\nif writer:\n    writer.close()\n# Kill the 7z pipe process.  terminate()+wait() can deadlock when the 7z process\n# is blocked on a full pipe write; kill() (SIGKILL) bypasses that.\nproc.kill()\nproc.wait()\nif SMOKE_MODE:\n    print(f\"train (SMOKE): {total_rows:,} rows  elapsed={_time.time()-t_train:.0f}s\")\nelse:\n    print(f\"train:        {total_rows:,} rows  elapsed={_time.time()-t_train:.0f}s\")\n\nprint(f\"\\nCSV\u2192parquet complete in {_time.time()-t0:.0f}s\")\n\nif SMOKE_MODE:\n    # Create data/dev/ slice \u2014 2 stores only \u2014 to cap pair universe for smoke test.\n    # Smoke cutoff (2013-01-15) predates test manifest dates (2017-08-16+), so\n    # preds won't join; submission fills 0 and row-count assert still runs.\n    _t_dev = _time.time()\n    _test_full  = pd.read_parquet(DATA / \"test.parquet\")\n    SMOKE_STORES = sorted(_test_full[\"store_nbr\"].unique())[:2]\n    print(f\"SMOKE: using stores {SMOKE_STORES} for dev pair-universe cap\")\n    os.makedirs(DATA / \"dev\", exist_ok=True)\n    _test_dev = _test_full[_test_full[\"store_nbr\"].isin(SMOKE_STORES)].copy()\n    _test_dev.to_parquet(DATA / \"dev/test.parquet\", index=False)\n    _train_full = pd.read_parquet(DATA / \"train.parquet\")   # only 2 chunks in SMOKE_MODE\n    _train_dev  = _train_full[_train_full[\"store_nbr\"].isin(SMOKE_STORES)].copy()\n    _train_dev.to_parquet(DATA / \"dev/train.parquet\", index=False)\n    print(\n        f\"SMOKE: dev slices written \"\n        f\"test={len(_test_dev):,}  train={len(_train_dev):,}  \"\n        f\"t={_time.time()-_t_dev:.1f}s\"\n    )"}, {"cell_type": "markdown", "metadata": {}, "source": "## Package setup\n\nThe next four cells write Python packages to `/kaggle/working/` using the `%%writefile` cell magic. Kaggle kernels start with an empty working directory; this pattern bundles the pipeline dependencies inline so the notebook is self-contained with no external installs required.\n\n- **`dataguard/`** \u2014 cutoff-enforcing data access layer. All reads go through `load()`, which raises if any requested column would leak data from after the forecast cutoff. Exceptions are logged explicitly (see the reader source).\n- **`models/lgbm_stage3.py`** \u2014 the full 72-feature LightGBM pipeline (~900 lines). Feature engineering, training loop, importance export, and inference."}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "%%writefile /kaggle/working/dataguard/__init__.py\n\"\"\"\ndataguard \u2014 the only sanctioned way for pipeline code to read data.\n\nLEAKAGE CONTRACT\n----------------\nEvery call returns rows with date STRICTLY BEFORE the supplied cutoff.\nThree documented exceptions exist; all require allow_future=True:\n\nEXCEPTION 1 \u2014 test table / onpromotion (competition-sanctioned):\n  Retailers know their own promotional schedules in advance.  The\n  competition supplies `onpromotion` for the 16-day prediction window\n  in test.csv.  Callers building promo features for the forecast horizon\n  may load this table without date filtering.\n\nEXCEPTION 2 \u2014 holidays table (calendar is public knowledge):\n  Government and company holiday calendars are published in advance.\n  A real forecaster knows next week's holidays today.  Loading the\n  holidays table past the cutoff is not leakage.\n\nEXCEPTION 3 \u2014 train table / onpromotion for forecast horizon:\n  Retailers know their own promo schedule for the upcoming window, the\n  same rationale as Exception 1.  However, train.csv holds historical\n  unit_sales which must NEVER be accessible past the cutoff.  When\n  allow_future=True is used on \"train\", dataguard enforces two guards:\n    \u2022 date filter runs FORWARD (date >= cutoff, not < cutoff)\n    \u2022 only columns [store_nbr, item_nbr, date, onpromotion] are returned\n  unit_sales and every other column are permanently withheld for future\n  rows.  Requesting them via the columns= parameter raises ValueError.\n\n  This exception closes the train/test asymmetry: the backtest now uses\n  train's future onpromotion for horizon features (same as the Kaggle\n  submission, which reads them from test.csv via Exception 1).\n\nAll three exceptions use the same flag: allow_future=True.\nIt is rejected (ValueError) for all other tables \u2014 oil prices and\ntransaction counts are genuinely unknowable in advance.\n\nEXCEPTION 4 \u2014 test.csv scoring manifest (assortment is known at forecast time):\n  test.csv's set of (store, item, date) triples is prediction-time knowledge.\n  Kaggle releases test.csv at the start of the forecast window; a real\n  retailer knows its own assortment and which products it will sell in\n  which stores.  Using the scoring manifest to define:\n    \u2022 the training population (dense pivot includes all test pairs)\n    \u2022 hierarchy peer sets (F6 group-mean features aggregate over test pairs)\n    \u2022 the inference scoring loop\n  is therefore not leakage \u2014 it mirrors information available to any\n  operational forecasting system at runtime.\n\n  This exception requires no allow_future flag and no code change in\n  dataguard; it is a documented rationale for existing behaviour.\n  The only column from test.csv that carries forward-looking signal is\n  onpromotion (covered by Exception 1).  Membership alone carries none.\n\"\"\"\n\nfrom dataguard.reader import load, load_stores, load_items  # noqa: F401\n"}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "%%writefile /kaggle/working/dataguard/reader.py\n\"\"\"\ndataguard.reader \u2014 cutoff-enforcing parquet access layer.\n\nAll public functions enforce: returned rows have date < cutoff.\n\nDOCUMENTED LEAKAGE EXCEPTIONS (see __init__.py for full rationale).\nThree tables may be loaded for future dates via allow_future=True:\n\n  1. \"test\"     \u2014 competition-provided onpromotion flags for the 16-day\n                  prediction window.  Retailers know their own promo\n                  schedules in advance; the competition legitimately\n                  supplies this column in test.csv.\n\n  2. \"holidays\" \u2014 a published government/company calendar.  Real\n                  forecasters know next week's public holidays today.\n                  Loading holidays past the cutoff is NOT leakage.\n\n  3. \"train\"    \u2014 onpromotion for the forecast horizon only.\n                  Same known-in-advance rationale as exceptions 1 and 2.\n                  Returns ONLY rows with date >= cutoff AND ONLY the four\n                  safe columns: store_nbr, item_nbr, date, onpromotion.\n                  unit_sales and all other columns are never exposed for\n                  future rows.  Requesting them via columns= raises ValueError.\n\n                  KNOWN LIMITATION: train's future onpromotion rows cover\n                  all (store, item) pairs that ever appeared in training \u2014\n                  including those with zero sales during the forecast window.\n                  test.csv only lists scoring pairs (non-zero rows in the\n                  competition's eyes).  Backtest promo coverage therefore\n                  exceeds submission promo coverage slightly; this is\n                  accepted variance and does not constitute leakage.\n\nNo other table accepts allow_future=True.  Oil prices and transaction\ncounts are genuinely unknowable in advance and remain strictly filtered.\n\nEXCEPTION 4 \u2014 scoring manifest (no code change required, see __init__.py):\n  test.csv's set of (store, item, date) triples is prediction-time knowledge.\n  Using the manifest to define the training population, the scoring loop, or\n  hierarchy peer sets is not leakage.  Membership carries no forward-looking\n  signal; only the onpromotion column does (covered by Exception 1).\n\"\"\"\n\nfrom __future__ import annotations\n\nimport os\nfrom pathlib import Path\n\nimport pandas as pd\n\n_DATA_ROOT = Path(os.environ.get(\"FAVORITA_DATA_ROOT\", Path(__file__).parent.parent / \"data\"))\n\n# Tables that carry a \"date\" column and must be cutoff-filtered.\n_DATED_TABLES = {\"train\", \"test\", \"oil\", \"transactions\", \"holidays\"}\n\n# Tables without a date column \u2014 returned as-is.\n_STATIC_TABLES = {\"stores\", \"items\"}\n\n# Tables for which allow_future=True bypasses the date filter entirely\n# (returns all rows regardless of date).\n_ALLOW_FUTURE_TABLES = {\"test\", \"holidays\"}\n\n# For \"train\" with allow_future=True: only these columns may be returned.\n# unit_sales and every other field are withheld \u2014 they represent unknowable\n# future ground truth.\n_TRAIN_PROMO_COLS: frozenset[str] = frozenset(\n    {\"store_nbr\", \"item_nbr\", \"date\", \"onpromotion\"}\n)\n\n\ndef _parquet_path(table: str, dev: bool) -> Path:\n    if dev:\n        path = _DATA_ROOT / \"dev\" / f\"{table}.parquet\"\n        if not path.exists():\n            # fall back to full table for static lookups and non-dev-sliced tables\n            path = _DATA_ROOT / f\"{table}.parquet\"\n    else:\n        path = _DATA_ROOT / f\"{table}.parquet\"\n    if not path.exists():\n        raise FileNotFoundError(f\"Parquet not found: {path}\")\n    return path\n\n\ndef load(\n    table: str,\n    cutoff: pd.Timestamp | str | None = None,\n    *,\n    dev: bool = True,\n    allow_future: bool = False,\n    columns: list[str] | None = None,\n) -> pd.DataFrame:\n    \"\"\"\n    Return a DataFrame for `table`.\n\n    Parameters\n    ----------\n    table : str\n        One of: train, test, oil, transactions, holidays, stores, items.\n    cutoff : Timestamp or str\n        Upper bound (exclusive) for normal access: rows with date >= cutoff\n        are dropped.  For allow_future=True on \"train\", used as LOWER bound\n        (inclusive): only rows with date >= cutoff are returned.\n        Required for all dated tables.\n    dev : bool\n        When True, reads from data/dev/ (2-store slice).  Set False only\n        for the final Kaggle notebook run.\n    allow_future : bool\n        LEAKAGE EXCEPTION \u2014 see __init__.py for full rationale.\n        Legal only for \"test\", \"holidays\", and \"train\" (restricted):\n\n          \"test\" / \"holidays\" : date filter is skipped entirely.\n          \"train\"             : date filter runs FORWARD (date >= cutoff);\n                               only _TRAIN_PROMO_COLS are returned.\n\n        Raises ValueError for any other table (oil, transactions).\n    columns : list[str] or None\n        Column whitelist.  Only active for allow_future=True on \"train\".\n        Every requested column must be in _TRAIN_PROMO_COLS; any column\n        outside that set (e.g. \"unit_sales\") raises ValueError.\n        Pass None (the default) to receive all four promo columns.\n    \"\"\"\n    # \u2500\u2500 validate allow_future \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    if allow_future and table not in _ALLOW_FUTURE_TABLES and table != \"train\":\n        raise ValueError(\n            f\"allow_future is only valid for tables \"\n            f\"{sorted(_ALLOW_FUTURE_TABLES | {'train'})}. \"\n            f\"Called with table='{table}'. \"\n            \"Oil and transactions are genuinely unknowable in advance.\"\n        )\n\n    # \u2500\u2500 validate columns= \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    # columns= is only meaningful for the train future-promo path.\n    if columns is not None:\n        if not (allow_future and table == \"train\"):\n            raise ValueError(\n                \"The columns= parameter is only valid for \"\n                \"allow_future=True on table='train'.\"\n            )\n        invalid = sorted(set(columns) - _TRAIN_PROMO_COLS)\n        if invalid:\n            raise ValueError(\n                f\"allow_future on 'train' exposes only \"\n                f\"{sorted(_TRAIN_PROMO_COLS)}. \"\n                f\"Columns {invalid} are not available for future dates \u2014 \"\n                \"they may contain target leakage (e.g. unit_sales).\"\n            )\n\n    # \u2500\u2500 read parquet \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    path = _parquet_path(table, dev=dev)\n    df = pd.read_parquet(path)\n\n    # \u2500\u2500 static tables \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    if table in _STATIC_TABLES:\n        return df\n\n    # \u2500\u2500 EXCEPTION 3: train future-promo \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    if allow_future and table == \"train\":\n        if cutoff is None:\n            raise ValueError(\"cutoff is required when using allow_future on 'train'\")\n        cutoff_ts = pd.Timestamp(cutoff)\n        _cols = sorted(_TRAIN_PROMO_COLS)  # deterministic column order\n        result = (\n            df.loc[df[\"date\"] >= cutoff_ts, _cols]\n            .reset_index(drop=True)\n        )\n        return result if columns is None else result[sorted(columns)]\n\n    # \u2500\u2500 EXCEPTIONS 1 & 2: test / holidays \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    if allow_future:\n        # LEAKAGE EXCEPTION: return without date restriction.\n        # \"test\"     -> competition onpromotion schedule for forecast window.\n        # \"holidays\" -> known-in-advance calendar.\n        return df\n\n    # \u2500\u2500 normal cutoff-filtered access \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    if cutoff is None:\n        raise ValueError(f\"cutoff is required for table='{table}'\")\n\n    cutoff_ts = pd.Timestamp(cutoff)\n\n    if \"date\" not in df.columns:\n        raise KeyError(f\"Table '{table}' has no 'date' column; cannot apply cutoff.\")\n\n    return df.loc[df[\"date\"] < cutoff_ts].reset_index(drop=True)\n\n\ndef load_stores(dev: bool = True) -> pd.DataFrame:\n    \"\"\"Convenience wrapper \u2014 stores has no date column.\"\"\"\n    return load(\"stores\", dev=dev)\n\n\ndef load_items(dev: bool = True) -> pd.DataFrame:\n    \"\"\"Convenience wrapper \u2014 items has no date column.\"\"\"\n    return load(\"items\", dev=dev)\n"}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "%%writefile /kaggle/working/models/__init__.py\n# models package\n"}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "%%writefile /kaggle/working/models/lgbm_stage3.py\n\"\"\"\n72-feature LightGBM pipeline for Corporaci\u00f3n Favorita sales forecasting.\n\nTrains 16 per-horizon models (h=1\u202616) on a 56-day rolling window; features\nspan rolling statistics, trends, promotions, calendar, hierarchy aggregates,\nsparsity indicators, oil/transactions, and family-level promo lift ratios.\nSee the README cell at the top of the notebook for scores and honesty notes.\n\"\"\"\n\nfrom __future__ import annotations\n\n# Memory management: matrices trimmed via .copy() (not slice-views) to release backing\n# arrays; gc.collect() + malloc_trim(0) returns freed pages to OS.\n# is_earthquake pruned \u2014 zero by construction for all 2017 training cutoffs (quake: 2016-04-16).\n# Feature importances written to env var FAVORITA_IMPORTANCE_PATH; set \"\" to skip.\n\nimport warnings\nimport numpy as np\nimport pandas as pd\nimport lightgbm as lgb\nimport psutil\nimport gc\nimport time\nimport os\nimport ctypes as _ct\n_libc = _ct.CDLL(\"libc.so.6\")   # for malloc_trim(0) \u2014 forces freed pages back to OS\n\nfrom dataguard import load, load_items, load_stores\n\n# \u2500\u2500 Family 8: External features \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\nFAMILY: str = \"F9-family-promo-lift\"\nSTAGE_DESCRIPTION: str = (\n    \"Stage 3-F9 family_promo_lift(on/off ratio from pre-cutoff history)\"\n)\n\nLAGS: list[int] = [1, 7, 14, 21, 28, 35, 42, 49, 56]\nHORIZON: int = 16\nTRAIN_DAYS: int = 56   # days of training history per forecast origin\n\n_ROLL_WINDOWS:         list[int] = [7, 14, 30, 60, 140]   # full stats\n_COND_WINDOWS:         list[int] = [7, 14, 30]             # promo-conditional means\n_PROMO_SHARE_WINDOWS:  list[int] = [7, 14, 30]             # F4: trailing promo rate\n\n# History pad: must satisfy both lag56+HORIZON-1=71 AND 140d rolling+TRAIN_DAYS=211.\n# Binding: 211.  n_dates = 267.\n_HISTORY_PAD: int = 211\n\n# Minimum column index for trimming feature matrices. The sales pivot has _HISTORY_PAD=211\n# columns before the first training date; the latest safe origin for h=16 is 211-(16-1)=196.\n# All 16 horizon models share the same trimmed slice; this is independent of TRAIN_DAYS.\nTRAIN_COL_START: int = 196\n\n# onpromotion three-state encoding\nPROMO_NA_CODE: float = 2.0   # \"untracked\" \u2014 pre-2014, distinct from False(0) / True(1)\n\nLGB_PARAMS: dict = {\n    \"objective\": \"regression\",\n    \"metric\": \"rmse\",\n    \"n_estimators\": 200,\n    \"learning_rate\": 0.05,\n    \"num_leaves\": 31,\n    \"min_child_samples\": 20,\n    \"feature_fraction\": 0.9,\n    \"bagging_fraction\": 0.8,\n    \"bagging_freq\": 5,\n    \"verbose\": -1,\n    \"n_jobs\": -1,\n    \"seed\": 42,  # overridden per blend member in the training cell\n}\n\n_ROLL_FEAT_NAMES: list[str] = (\n    [f\"rm{w}\"       for w in _ROLL_WINDOWS]\n  + [f\"rmed{w}\"     for w in _ROLL_WINDOWS]\n  + [f\"rstd{w}\"     for w in _ROLL_WINDOWS]\n  + [f\"rmin{w}\"     for w in _ROLL_WINDOWS]\n  + [f\"rmax{w}\"     for w in _ROLL_WINDOWS]\n  + [f\"promo_mean{w}\"   for w in _COND_WINDOWS]\n  + [f\"nopromo_mean{w}\" for w in _COND_WINDOWS]\n)\n_ROLL_FEAT_NAMES = [f for f in _ROLL_FEAT_NAMES if f not in {\"rmin60\", \"rmin140\"}]  # rmin60/rmin140: bottom-decile gain; rmin7/14/30 cover the diminishing-window signal\n\n# Trend diffs, ratios, and slopes \u2014 all derived at the same pre-origin window.\n_TREND_PAIRS: list[tuple[int, int]] = [(7, 14), (7, 30), (7, 60), (14, 60), (30, 140)]\n_TREND_FEAT_NAMES: list[str] = (\n    [f\"d{a}_{b}\" for a, b in _TREND_PAIRS]\n  + [f\"r{a}_{b}\" for a, b in _TREND_PAIRS]\n  + [\"slope7\", \"slope14\"]\n)\n_TREND_FEAT_NAMES = [f for f in _TREND_FEAT_NAMES if f not in {\"d7_14\"}]  # d7_14: bottom-decile; d7_30/d7_60 cover the short-window trend\n\n_F4_FEAT_NAMES: list[str] = (\n    [f\"promo_share{w}\" for w in _PROMO_SHARE_WINDOWS]\n    + [\"days_since_promo\"]\n    + [\"promo_count_fwd16\"]\n)\n\n_F5_FEAT_NAMES: list[str] = [\n    \"dom\", \"days_to_payday\", \"days_since_nat_hol\",\n    # month/is_nat_holiday/days_to_nat_hol/is_local_holiday: bottom-decile or zero gain\n]\n\n_F6_FEAT_NAMES: list[str] = [\n    \"store_agg_rm30\", \"item_agg_rm30\",\n    \"family_agg_rm30\", \"class_agg_rm30\",\n    # cluster_agg_rm30: degenerate when 1 store/cluster (dev); family/class/store levels healthy\n]\n\n_F7_ZERO_WINDOWS: list[int] = [30, 140]\n_F7_FEAT_NAMES: list[str] = (\n    [\"days_since_first_sale\", \"days_since_last_sale\"]\n    + [f\"zero_share{w}\" for w in _F7_ZERO_WINDOWS]\n)\n\n_F8_FEAT_NAMES: list[str] = [\n    \"oil_rm30\", \"oil_r7_30\",\n    \"txn_rm30\", \"txn_r7_30\",\n]\n\n_F9_FEAT_NAMES: list[str] = [\"family_promo_lift\"]\n\nFEATURE_NAMES: list[str] = (\n    [f\"lag{k}\" for k in LAGS]\n    + [\"dow\", \"onpromotion\"]\n    + _ROLL_FEAT_NAMES\n    + _TREND_FEAT_NAMES\n    + _F4_FEAT_NAMES\n    + _F5_FEAT_NAMES\n    + _F6_FEAT_NAMES\n    + _F7_FEAT_NAMES\n    + _F8_FEAT_NAMES\n    + _F9_FEAT_NAMES\n)\n\n\ndef _rolling_train_mats(\n    sales: np.ndarray, promo: np.ndarray\n) -> dict[str, np.ndarray]:\n    \"\"\"\n    Precompute rolling stat matrices for training.\n\n    Returns dict: name \u2192 (n_pairs, n_dates) float32 where value at [p, j]\n    is the statistic of sales[p, j-w : j] (window ending strictly BEFORE j).\n    This matches the information available at forecast origin j.\n    \"\"\"\n    n_pairs, n_dates = sales.shape\n    # Transpose: pandas rolling operates along rows (axis=0).\n    s = pd.DataFrame(sales.T)  # (n_dates, n_pairs)\n\n    promo_mask = (promo == 1.0).astype(\"float32\")  # PROMO_NA_CODE=2 \u2192 not promo\n    pm = pd.DataFrame(promo_mask.T)\n    ps = pd.DataFrame((sales * promo_mask).T)\n    nm = pd.DataFrame((1.0 - promo_mask).T)        # 1 where promo=0 (no-promo day)\n    ns = pd.DataFrame((sales * (1.0 - promo_mask)).T)\n\n    out: dict[str, np.ndarray] = {}\n\n    for w in _ROLL_WINDOWS:\n        r = s.rolling(w, min_periods=1)\n        # .shift(1): shift results forward so position j holds the stat of [j-w, j)\n        out[f\"rm{w}\"]   = r.mean().shift(1).fillna(0.0).values.T.astype(\"float32\")\n        out[f\"rmed{w}\"] = r.median().shift(1).fillna(0.0).values.T.astype(\"float32\")\n        std = r.std(ddof=1).shift(1).fillna(0.0)\n        out[f\"rstd{w}\"] = std.values.T.astype(\"float32\")\n        out[f\"rmin{w}\"] = r.min().shift(1).fillna(0.0).values.T.astype(\"float32\")\n        out[f\"rmax{w}\"] = r.max().shift(1).fillna(0.0).values.T.astype(\"float32\")\n\n    for w in _COND_WINDOWS:\n        p_sum = ps.rolling(w, min_periods=1).sum().shift(1).fillna(0.0).values.T.astype(\"float32\")\n        p_cnt = pm.rolling(w, min_periods=1).sum().shift(1).fillna(0.0).values.T.astype(\"float32\")\n        out[f\"promo_mean{w}\"] = np.where(p_cnt > 0, p_sum / (p_cnt + 1e-9), 0.0).astype(\"float32\")\n\n        n_sum = ns.rolling(w, min_periods=1).sum().shift(1).fillna(0.0).values.T.astype(\"float32\")\n        n_cnt = nm.rolling(w, min_periods=1).sum().shift(1).fillna(0.0).values.T.astype(\"float32\")\n        out[f\"nopromo_mean{w}\"] = np.where(n_cnt > 0, n_sum / (n_cnt + 1e-9), 0.0).astype(\"float32\")\n\n    # \u2500\u2500 Trend: diffs, ratios, slopes \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    _eps = 1e-9\n\n    for a, b in _TREND_PAIRS:\n        rm_a = out[f\"rm{a}\"]\n        rm_b = out[f\"rm{b}\"]\n        out[f\"d{a}_{b}\"] = (rm_a - rm_b).astype(\"float32\")\n        out[f\"r{a}_{b}\"] = np.clip(rm_a / (rm_b + _eps), 0.0, 5.0).astype(\"float32\")\n\n    for sw in [7, 14]:\n        x_s = np.arange(sw, dtype=\"float32\") - (sw - 1) / 2.0   # centred\n        x2  = float((x_s ** 2).sum())\n        slope_m = np.zeros((n_pairs, n_dates), dtype=\"float32\")\n        for j in range(sw, n_dates):\n            win = sales[:, j - sw : j]   # (n_pairs, sw) \u2014 strictly before j\n            slope_m[:, j] = (win * x_s).sum(axis=1) / x2\n        out[f\"slope{sw}\"] = slope_m\n\n    # \u2500\u2500 F4: Promo share \u2014 fraction of trailing w days on promotion (pre-origin) \u2500\n    for w in _PROMO_SHARE_WINDOWS:\n        out[f\"promo_share{w}\"] = (\n            pm.rolling(w, min_periods=1).mean().shift(1).fillna(0.0)\n            .values.T.astype(\"float32\")\n        )\n\n    # \u2500\u2500 F4: Days since last promo (pre-origin) \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    # Value at [p, j] = j minus last position k < j where promo==1, capped at TRAIN_DAYS.\n    row_ix = np.where(\n        pm.values == 1.0,\n        np.arange(n_dates, dtype=\"float64\").reshape(-1, 1),  # (n_dates, 1) broadcasts\n        np.nan,\n    )\n    last_promo = pd.DataFrame(row_ix).ffill().shift(1)  # strictly before j\n    curr = np.arange(n_dates, dtype=\"float64\").reshape(-1, 1)\n    dsp = curr - last_promo.values\n    dsp = np.where(np.isnan(dsp), float(TRAIN_DAYS), dsp)\n    out[\"days_since_promo\"] = np.clip(dsp, 0.0, float(TRAIN_DAYS)).T.astype(\"float32\")\n\n    # \u2500\u2500 F7: Sparsity \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    nonzero = (sales > 0)   # (n_pairs, n_dates) bool\n    j_idx   = np.arange(n_dates, dtype=\"float64\")\n\n    # days_since_first_sale[p, j]: strictly before j\n    has_any   = nonzero.any(axis=1)   # (n_pairs,)\n    first_pos = np.where(has_any, np.argmax(nonzero, axis=1), n_dates).astype(\"float64\")\n    dsf = np.where(\n        first_pos[:, None] < j_idx[None, :],\n        j_idx[None, :] - first_pos[:, None],\n        j_idx[None, :],\n    )\n    out[\"days_since_first_sale\"] = np.clip(dsf, 0, n_dates).astype(\"float32\")\n\n    # days_since_last_sale[p, j]: last nonzero strictly before j, cap TRAIN_DAYS\n    last_pos_mat = (\n        pd.DataFrame(np.where(nonzero.T, j_idx[:, None], np.nan))\n        .ffill().shift(1).values.T   # (n_pairs, n_dates)\n    )\n    dsl = np.where(\n        np.isnan(last_pos_mat),\n        float(TRAIN_DAYS),\n        j_idx[None, :] - last_pos_mat,\n    )\n    out[\"days_since_last_sale\"] = np.clip(dsl, 0, TRAIN_DAYS).astype(\"float32\")\n\n    # zero_share: fraction of trailing w days with zero sales (pre-origin)\n    zi = pd.DataFrame((sales == 0).astype(\"float32\").T)   # (n_dates, n_pairs)\n    for w in _F7_ZERO_WINDOWS:\n        out[f\"zero_share{w}\"] = (\n            zi.rolling(w, min_periods=1).mean().shift(1).fillna(1.0)\n            .values.T.astype(\"float32\")\n        )\n\n    return out\n\n\ndef _rolling_inf_vecs(\n    sales: np.ndarray, promo: np.ndarray, n_dates: int\n) -> dict[str, np.ndarray]:\n    \"\"\"\n    Compute rolling stats at the forecast origin (= cutoff) for inference.\n    Uses the last w columns of the pivot (= last w days before cutoff).\n    Returns dict: name \u2192 (n_pairs,) float32.\n    \"\"\"\n    promo_mask = (promo == 1.0).astype(\"float32\")\n    out: dict[str, np.ndarray] = {}\n\n    for w in _ROLL_WINDOWS:\n        sl = sales[:, max(0, n_dates - w) : n_dates]  # (n_pairs, \u2264w)\n        out[f\"rm{w}\"]   = sl.mean(axis=1).astype(\"float32\")\n        out[f\"rmed{w}\"] = np.median(sl, axis=1).astype(\"float32\")\n        out[f\"rstd{w}\"] = (\n            np.nan_to_num(sl.std(axis=1, ddof=1), nan=0.0).astype(\"float32\")\n            if sl.shape[1] > 1 else np.zeros(sl.shape[0], \"float32\")\n        )\n        out[f\"rmin{w}\"] = sl.min(axis=1).astype(\"float32\")\n        out[f\"rmax{w}\"] = sl.max(axis=1).astype(\"float32\")\n\n    for w in _COND_WINDOWS:\n        sl  = sales[:, max(0, n_dates - w) : n_dates]\n        pm  = promo_mask[:, max(0, n_dates - w) : n_dates]\n        nm  = (1.0 - pm)\n        p_cnt = pm.sum(axis=1)\n        n_cnt = nm.sum(axis=1)\n        out[f\"promo_mean{w}\"]   = np.where(p_cnt > 0, (sl * pm).sum(axis=1) / (p_cnt + 1e-9), 0.0).astype(\"float32\")\n        out[f\"nopromo_mean{w}\"] = np.where(n_cnt > 0, (sl * nm).sum(axis=1) / (n_cnt + 1e-9), 0.0).astype(\"float32\")\n\n    # \u2500\u2500 Trend: diffs, ratios, slopes \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    _eps = 1e-9\n    for a, b in _TREND_PAIRS:\n        rm_a = out[f\"rm{a}\"]\n        rm_b = out[f\"rm{b}\"]\n        out[f\"d{a}_{b}\"] = (rm_a - rm_b).astype(\"float32\")\n        out[f\"r{a}_{b}\"] = np.clip(rm_a / (rm_b + _eps), 0.0, 5.0).astype(\"float32\")\n\n    for sw in [7, 14]:\n        x_s = np.arange(sw, dtype=\"float32\") - (sw - 1) / 2.0\n        x2  = float((x_s ** 2).sum())\n        sl  = sales[:, max(0, n_dates - sw) : n_dates]   # (n_pairs, sw)\n        out[f\"slope{sw}\"] = ((sl * x_s).sum(axis=1) / x2).astype(\"float32\")\n\n    # \u2500\u2500 F4: Promo share (inference) \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    for w in _PROMO_SHARE_WINDOWS:\n        pm_w = promo_mask[:, max(0, n_dates - w) : n_dates]\n        out[f\"promo_share{w}\"] = pm_w.mean(axis=1).astype(\"float32\")\n\n    # \u2500\u2500 F4: Days since last promo (inference) \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    # n_dates columns \u2192 last column is \"yesterday\" relative to cutoff.\n    # days_since = (n_dates - 1 - last_promo_pos) + 1 = n_dates - last_promo_pos.\n    has_promo = promo_mask.any(axis=1)                         # (n_pairs,) bool\n    k_from_right = np.argmax(promo_mask[:, ::-1], axis=1)     # 0 \u2192 yesterday was promo\n    dsp_inf = np.where(has_promo, 1 + k_from_right, TRAIN_DAYS)\n    out[\"days_since_promo\"] = np.clip(dsp_inf, 0, TRAIN_DAYS).astype(\"float32\")\n    # NOTE: promo_count_fwd16 requires horizon_dates; computed in lgbm_stage3().\n\n    # \u2500\u2500 F7: Sparsity (inference) \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    nonzero_inf = (sales > 0)   # (n_pairs, n_dates)\n    has_any_inf = nonzero_inf.any(axis=1)\n\n    first_pos_inf = np.where(has_any_inf, np.argmax(nonzero_inf, axis=1), n_dates)\n    out[\"days_since_first_sale\"] = np.clip(\n        np.where(has_any_inf, n_dates - first_pos_inf, n_dates).astype(\"float32\"),\n        0, n_dates,\n    )\n\n    # k_from_right_sale: 0 \u2192 yesterday was a nonzero sale\n    k_from_right_sale = np.argmax(nonzero_inf[:, ::-1], axis=1)\n    out[\"days_since_last_sale\"] = np.clip(\n        np.where(has_any_inf, 1 + k_from_right_sale, TRAIN_DAYS).astype(\"float32\"),\n        0, TRAIN_DAYS,\n    )\n\n    for w in _F7_ZERO_WINDOWS:\n        sl = sales[:, max(0, n_dates - w) : n_dates]\n        out[f\"zero_share{w}\"] = (sl == 0).mean(axis=1).astype(\"float32\")\n\n    return out\n\n\n# \u2500\u2500 F5: Calendar helpers \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n_HOL_TYPES    = frozenset({\"Holiday\", \"Transfer\", \"Additional\", \"Bridge\"})\n_QUAKE_START  = pd.Timestamp(\"2016-04-16\").date()\n_QUAKE_END    = pd.Timestamp(\"2016-05-16\").date()\n_HOL_DIST_CAP = 5   # max days reported to/from nearest national holiday\n\n\ndef _load_holiday_sets(hol: pd.DataFrame) -> tuple[set, set, set]:\n    \"\"\"\n    Parse holidays DataFrame into three frozensets used by calendar features.\n    Returns: (nat_set, reg_set, loc_set).\n      nat_set: {date} \u2014 national holiday observances\n      reg_set: {(locale_name, date)} \u2014 regional holiday observances\n      loc_set: {(locale_name, date)} \u2014 local holiday observances\n    Earthquake Events are excluded (flagged via _QUAKE_START/_QUAKE_END instead).\n    \"\"\"\n    obs = hol[hol[\"type\"].isin(_HOL_TYPES)]\n    nat_set = set(obs[obs[\"locale\"] == \"National\"][\"date\"].dt.date)\n    reg = obs[obs[\"locale\"] == \"Regional\"]\n    reg_set = set(zip(reg[\"locale_name\"].astype(str), reg[\"date\"].dt.date))\n    loc = obs[obs[\"locale\"] == \"Local\"]\n    loc_set = set(zip(loc[\"locale_name\"].astype(str), loc[\"date\"].dt.date))\n    return nat_set, reg_set, loc_set\n\n\ndef _nat_hol_dists(\n    dates: pd.DatetimeIndex, nat_set: set, cap: int = _HOL_DIST_CAP\n) -> tuple[np.ndarray, np.ndarray]:\n    \"\"\"\n    For each date return (days_to_next_holiday, days_since_last_holiday).\n    Both capped at `cap`.  Value = cap when no holiday found within cap days.\n    Uses ordinal arithmetic for speed; O(n \u00d7 log H) via searchsorted.\n    \"\"\"\n    nat_ord = np.array(sorted(d.toordinal() for d in nat_set))\n    d_ord   = np.array([d.date().toordinal() for d in dates])\n    n = len(d_ord)\n\n    idx = np.searchsorted(nat_ord, d_ord)  # insertion point in sorted nat_ord\n\n    # days to next: nat_ord[idx] - d_ord (idx < len \u2192 holiday exists at or after d)\n    days_to = np.full(n, float(cap), dtype=\"float32\")\n    has_next = idx < len(nat_ord)\n    days_to[has_next] = np.minimum(\n        (nat_ord[idx[has_next]] - d_ord[has_next]).astype(\"float32\"), cap\n    )\n\n    # days since last: d_ord - nat_ord[idx-1] (idx > 0 \u2192 holiday exists at or before d)\n    days_from = np.full(n, float(cap), dtype=\"float32\")\n    has_prev = idx > 0\n    days_from[has_prev] = np.minimum(\n        (d_ord[has_prev] - nat_ord[idx[has_prev] - 1]).astype(\"float32\"), cap\n    )\n\n    return days_to, days_from\n\n\ndef _payday_prox(dates: pd.DatetimeIndex, cap: int = 7) -> np.ndarray:\n    \"\"\"Min(|dom\u22121|, |dom\u221215|, last_day\u2212dom) for each date, capped at cap.\"\"\"\n    out = np.empty(len(dates), dtype=\"float32\")\n    for i, d in enumerate(dates):\n        dom  = d.day\n        last = d.days_in_month\n        out[i] = float(min(dom - 1, abs(dom - 15), last - dom, cap))\n    return out\n\n\ndef _hier_agg_mat(mat: np.ndarray, groups: np.ndarray) -> np.ndarray:\n    \"\"\"\n    (n_pairs, n_dates): value at [p, j] = mean(mat[peers, j]) where peers share groups[p].\n    Singleton groups get their own value (no information leakage \u2014 the value equals mat[p]).\n    \"\"\"\n    out = np.empty_like(mat)\n    for g in pd.unique(groups):\n        mask = groups == g\n        out[mask] = mat[mask].mean(axis=0, keepdims=True)\n    return out\n\n\ndef _hier_agg_vec(vec: np.ndarray, groups: np.ndarray) -> np.ndarray:\n    \"\"\"(n_pairs,): value at [p] = mean(vec[peers]) where peers share groups[p].\"\"\"\n    out = np.empty_like(vec)\n    for g in pd.unique(groups):\n        mask = groups == g\n        out[mask] = vec[mask].mean()\n    return out\n\n\ndef _promo_float(col: pd.Series) -> np.ndarray:\n    \"\"\"Encode boolean/NA onpromotion: True\u21921.0, False\u21920.0, NA\u2192PROMO_NA_CODE.\"\"\"\n    not_na = col.notna().values\n    vals   = col.fillna(False).astype(bool).values\n    result = np.full(len(col), PROMO_NA_CODE, dtype=\"float32\")\n    result[not_na & vals]  = 1.0\n    result[not_na & ~vals] = 0.0\n    return result\n\n\ndef lgbm_stage3(\n    cutoff: pd.Timestamp,\n    horizon_dates: pd.DatetimeIndex,\n    dev: bool = True,\n) -> pd.DataFrame:\n    \"\"\"\n    Train 16 per-lead-time LightGBM models and return predictions.\n\n    Compatible with eval.backtest.backtest() signature:\n      returns DataFrame[store_nbr, item_nbr, date, unit_sales_pred].\n    \"\"\"\n    _t0 = time.time()   # wall-clock anchor for phase prints\n    cutoff = pd.Timestamp(cutoff)\n\n    # \u2500\u2500 1. load data \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    train_df = load(\"train\", cutoff=cutoff, dev=dev)\n    items_df  = load_items(dev=dev)\n    # EXCEPTION 1 (dataguard v1.0): test.csv onpromotion \u2014 covers submission window.\n    test_df   = load(\"test\", dev=dev, allow_future=True)\n    # EXCEPTION 3 (dataguard v1.1): train future promo for backtest horizons.\n    # Empty when cutoff >= train end date (submission path).\n    train_future_df = load(\"train\", cutoff=cutoff, dev=dev, allow_future=True)\n\n    perishable_items: set = set(\n        items_df.loc[items_df[\"perishable\"] == 1, \"item_nbr\"]\n    )\n\n    scoring_pairs = (\n        test_df[[\"store_nbr\", \"item_nbr\"]]\n        .drop_duplicates()\n        .reset_index(drop=True)\n    )\n    n_pairs = len(scoring_pairs)\n    pair_mi = pd.MultiIndex.from_frame(scoring_pairs)\n    print(f\"[phase 0] data loaded: n_pairs={n_pairs:,}  cutoff={cutoff.date()}  \"\n          f\"t={time.time()-_t0:.1f}s\", flush=True)\n\n    # \u2500\u2500 2. build pivot matrices \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    history_days = TRAIN_DAYS + _HISTORY_PAD   # 127\n    pivot_start  = cutoff - pd.Timedelta(days=history_days)\n    pivot_dates  = pd.date_range(pivot_start, cutoff - pd.Timedelta(days=1), freq=\"D\")\n    n_dates      = len(pivot_dates)   # 127\n\n    train_df = train_df.copy()\n    train_df[\"unit_sales\"] = train_df[\"unit_sales\"].clip(lower=0.0)\n    train_df[\"onpromo_f\"]  = _promo_float(train_df[\"onpromotion\"])\n\n    train_win = train_df[train_df[\"date\"] >= pivot_start].merge(\n        scoring_pairs, on=[\"store_nbr\", \"item_nbr\"]\n    )\n\n    def _pivot(col: str, default: float = 0.0) -> np.ndarray:\n        \"\"\"(n_pairs, n_dates) float32; missing cells \u2192 default.\"\"\"\n        agg = \"sum\" if col == \"unit_sales\" else \"max\"\n        wide = (\n            train_win.pivot_table(\n                index=[\"store_nbr\", \"item_nbr\"],\n                columns=\"date\",\n                values=col,\n                aggfunc=agg,\n            )\n            .reindex(index=pair_mi, columns=pivot_dates, fill_value=default)\n            .fillna(default)\n            .astype(\"float32\")\n        )\n        return wide.values\n\n    sales_mat = _pivot(\"unit_sales\", default=0.0)\n    # default=0.0: pairs absent from train_win on a date are treated as no-promo,\n    # NOT as PROMO_NA_CODE.  PROMO_NA_CODE is reserved for actual rows in train.csv\n    # where onpromotion IS NULL (pre-2014 untracked regime).  Using PROMO_NA_CODE\n    # as the absent-pair default would leak future transaction status: the model\n    # would learn \"PROMO_NA_CODE \u2192 y=0\" from training, then at inference the\n    # absence of a pair in train_future_df would be an illicit zero-sales signal.\n    promo_mat = _pivot(\"onpromo_f\",  default=0.0)\n    print(f\"[phase 1] pivot built: sales_mat {sales_mat.shape}  \"\n          f\"t={time.time()-_t0:.1f}s\", flush=True)\n\n    pair_weights = np.array(\n        [1.25 if item in perishable_items else 1.0\n         for _, item in scoring_pairs.itertuples(index=False)],\n        dtype=\"float32\",\n    )\n\n    # \u2500\u2500 2b. precompute rolling stat matrices \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    roll_train = _rolling_train_mats(sales_mat, promo_mat)\n    # RAM: .copy() not slice-view \u2014 views keep backing data pinned in RAM\n    roll_train = {k: v[:, TRAIN_COL_START:].copy() for k, v in roll_train.items()}\n    gc.collect(); _libc.malloc_trim(0)   # return freed pages to OS\n    _nc = next(iter(roll_train.values())).shape[1]\n    print(f\"[phase 2] roll_train: {len(roll_train)} mats \u00d7 {_nc} cols  \"\n          f\"RAM={psutil.Process().memory_info().rss/1e9:.2f}GB  t={time.time()-_t0:.1f}s\", flush=True)\n\n    # \u2500\u2500 2c. F5: Calendar precomputation \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    hol_df    = load(\"holidays\", dev=dev, allow_future=True)\n    stores_df = load_stores(dev=dev)\n    nat_set, reg_set, loc_set = _load_holiday_sets(hol_df)\n    store_city  = stores_df.set_index(\"store_nbr\")[\"city\"].to_dict()\n    store_state = stores_df.set_index(\"store_nbr\")[\"state\"].to_dict()\n\n    # Scalar calendar vectors over pivot_dates (n_dates,)\n    pivot_dom    = pivot_dates.day.values.astype(\"float32\")\n    pivot_month  = pivot_dates.month.values.astype(\"float32\")\n    pivot_payday = _payday_prox(pivot_dates)\n    pivot_nat_hol = np.array([1.0 if d.date() in nat_set else 0.0 for d in pivot_dates], \"float32\")\n    pivot_quake   = np.array(\n        [1.0 if _QUAKE_START <= d.date() <= _QUAKE_END else 0.0 for d in pivot_dates], \"float32\"\n    )\n    pivot_days_to, pivot_days_from = _nat_hol_dists(pivot_dates, nat_set)\n\n    # Per-pair store holiday arrays (n_pairs, n_dates) \u2014 identical within same store\n    _unique_stores = scoring_pairs[\"store_nbr\"].unique()\n    _store_reg = {\n        s: np.array([1.0 if (store_state[s], d.date()) in reg_set else 0.0 for d in pivot_dates], \"float32\")\n        for s in _unique_stores\n    }\n    _store_loc = {\n        s: np.array([1.0 if (store_city[s], d.date()) in loc_set else 0.0 for d in pivot_dates], \"float32\")\n        for s in _unique_stores\n    }\n    pair_stores_arr = scoring_pairs[\"store_nbr\"].values\n    pair_reg_hol = np.vstack([_store_reg[s] for s in pair_stores_arr])   # (n_pairs, n_dates)\n    pair_loc_hol = np.vstack([_store_loc[s] for s in pair_stores_arr])   # (n_pairs, n_dates)\n    print(f\"[phase 3] F5 calendar  RAM={psutil.Process().memory_info().rss/1e9:.2f}GB  t={time.time()-_t0:.1f}s\", flush=True)\n\n    # \u2500\u2500 2d. F6: Hierarchy aggregate matrices \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    # Group means of rm30 at each hierarchy level; same pre-origin guarantee as rm30.\n    grp_store   = scoring_pairs[\"store_nbr\"].values\n    grp_item    = scoring_pairs[\"item_nbr\"].values\n    grp_family  = scoring_pairs[\"item_nbr\"].map(\n        items_df.set_index(\"item_nbr\")[\"family\"]\n    ).values\n    grp_class   = scoring_pairs[\"item_nbr\"].map(\n        items_df.set_index(\"item_nbr\")[\"class\"]\n    ).values\n    grp_cluster = scoring_pairs[\"store_nbr\"].map(\n        stores_df.set_index(\"store_nbr\")[\"cluster\"]\n    ).values\n    rm30_mat = roll_train[\"rm30\"]\n    hier_mats = {\n        \"store_agg_rm30\":   _hier_agg_mat(rm30_mat, grp_store),\n        \"item_agg_rm30\":    _hier_agg_mat(rm30_mat, grp_item),\n        \"family_agg_rm30\":  _hier_agg_mat(rm30_mat, grp_family),\n        \"class_agg_rm30\":   _hier_agg_mat(rm30_mat, grp_class),\n        \"cluster_agg_rm30\": _hier_agg_mat(rm30_mat, grp_cluster),\n    }\n    print(f\"[phase 4] F6 hierarchy  RAM={psutil.Process().memory_info().rss/1e9:.2f}GB  t={time.time()-_t0:.1f}s\", flush=True)\n\n    # \u2500\u2500 2g. F8: External feature matrices (oil + transactions) \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    # Oil: global scalar series; NaN (weekends/holidays) forward-filled then backward-filled.\n    oil_raw = (\n        load(\"oil\", cutoff=cutoff, dev=dev)\n        .set_index(\"date\")[\"dcoilwtico\"]\n        .reindex(pivot_dates)\n        .ffill().bfill()\n        .fillna(50.0)   # fallback \u2014 shouldn't be reached in 2017 window\n        .values.astype(\"float32\")\n    )\n    _oil_s       = pd.Series(oil_raw)\n    # shift(1) makes rm30[t] = trailing mean ending at t-1 (strictly pre-origin).\n    # ffill() rather than bfill() so the boundary NaN at position 0 is filled\n    # from before (earliest known value), not from position 1 (a future value).\n    _oil_rm30_v  = _oil_s.rolling(30, min_periods=1).mean().shift(1).ffill().fillna(oil_raw[0]).values.astype(\"float32\")\n    _oil_rm7_v   = _oil_s.rolling(7,  min_periods=1).mean().shift(1).ffill().fillna(oil_raw[0]).values.astype(\"float32\")\n    _oil_r7_30_v = np.clip(_oil_rm7_v / (_oil_rm30_v + 1e-9), 0.0, 5.0).astype(\"float32\")\n    f8_oil_rm30_mat  = np.tile(_oil_rm30_v,  (n_pairs, 1))\n    f8_oil_r7_30_mat = np.tile(_oil_r7_30_v, (n_pairs, 1))\n\n    # Transactions: per-store daily counts; missing dates forward-filled then zero-filled.\n    txn_raw = load(\"transactions\", cutoff=cutoff, dev=dev)\n    _ust = scoring_pairs[\"store_nbr\"].unique()   # unique stores in dev slice\n    txn_wide = (\n        txn_raw.pivot_table(index=\"store_nbr\", columns=\"date\", values=\"transactions\", aggfunc=\"sum\")\n        .reindex(index=_ust, columns=pivot_dates)\n        .ffill(axis=1).fillna(0.0)\n        .astype(\"float32\")\n    )\n    _ust_pos = {s: i for i, s in enumerate(_ust)}\n    _pair_spos = np.array([_ust_pos[s] for s in scoring_pairs[\"store_nbr\"].values])\n    _txn_rm30_store = np.zeros((len(_ust), n_dates), dtype=\"float32\")\n    _txn_rm7_store  = np.zeros_like(_txn_rm30_store)\n    for _si in range(len(_ust)):\n        _ts = pd.Series(txn_wide.iloc[_si].values)\n        _txn_rm30_store[_si] = _ts.rolling(30, min_periods=1).mean().shift(1).fillna(0.0).values\n        _txn_rm7_store[_si]  = _ts.rolling(7,  min_periods=1).mean().shift(1).fillna(0.0).values\n    f8_txn_rm30_mat  = _txn_rm30_store[_pair_spos]   # (n_pairs, n_dates)\n    f8_txn_r7_30_mat = np.clip(\n        _txn_rm7_store[_pair_spos] / (_txn_rm30_store[_pair_spos] + 1e-9), 0.0, 5.0\n    ).astype(\"float32\")\n\n    f8_mats: dict[str, np.ndarray] = {\n        \"oil_rm30\":  f8_oil_rm30_mat,\n        \"oil_r7_30\": f8_oil_r7_30_mat,\n        \"txn_rm30\":  f8_txn_rm30_mat,\n        \"txn_r7_30\": f8_txn_r7_30_mat,\n    }\n    # RAM: .copy() releases 4 \u00d7 225 MB backing; gc + malloc_trim return pages to OS\n    f8_mats = {k: v[:, TRAIN_COL_START:].copy() for k, v in f8_mats.items()}\n    gc.collect(); _libc.malloc_trim(0)\n    print(f\"[phase 5] F8 oil+txn  RAM={psutil.Process().memory_info().rss/1e9:.2f}GB  t={time.time()-_t0:.1f}s\", flush=True)\n\n    # \u2500\u2500 2h. F9: Family promo lift \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    # Ratio of mean unit_sales on-promo vs off-promo from pre-cutoff history.\n    # Excludes pre-2014 rows (onpromotion untracked = NaN \u2192 onpromo_f = 2.0).\n    _lift_source = train_df[train_df[\"onpromotion\"].notna()].merge(\n        items_df[[\"item_nbr\", \"family\"]], on=\"item_nbr\", how=\"left\"\n    )\n    _lift_grp = (\n        _lift_source.groupby([\"family\", \"onpromotion\"])[\"unit_sales\"]\n        .mean()\n        .unstack(fill_value=np.nan)\n    )\n    _family_lift_map: dict = {}\n    for _fam in _lift_grp.index:\n        _off = _lift_grp.loc[_fam].get(False, np.nan)\n        _on  = _lift_grp.loc[_fam].get(True,  np.nan)\n        if np.isnan(float(_off)) or float(_off) < 0.01 or np.isnan(float(_on)):\n            _family_lift_map[_fam] = 1.0\n        else:\n            _family_lift_map[_fam] = float(np.clip(float(_on) / float(_off), 0.1, 10.0))\n    _pair_families = scoring_pairs[\"item_nbr\"].map(\n        items_df.set_index(\"item_nbr\")[\"family\"]\n    ).values\n    f9_lift_vec = np.array(\n        [_family_lift_map.get(f, 1.0) for f in _pair_families], dtype=\"float32\"\n    )\n    f9_lift_mat = np.tile(f9_lift_vec[:, np.newaxis], (1, n_dates))   # (n_pairs, n_dates)\n    # RAM: .copy() releases 225 MB tile backing\n    f9_lift_mat = f9_lift_mat[:, TRAIN_COL_START:].copy()\n    print(f\"[phase 6] F9 promo_lift  RAM={psutil.Process().memory_info().rss/1e9:.2f}GB  t={time.time()-_t0:.1f}s\", flush=True)\n\n    # \u2500\u2500 2e. precompute forward promo count (F4 ALLOWED EXCEPTION) \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    # promo_fwd16[p, j] = count of promo days in promo_mat[p, j : j+HORIZON].\n    # For training origins within the pivot, uses actual promo_mat values.\n    # For origins near pivot end, window is truncated (defaults to 0 for missing days).\n    # This is the competition-provided future onpromotion \u2014 never unit_sales.\n    promo_cs = np.cumsum(promo_mat, axis=1)   # (n_pairs, n_dates) prefix sums\n    end_idx  = np.minimum(np.arange(n_dates) + HORIZON, n_dates) - 1  # (n_dates,)\n    promo_fwd16_mat = promo_cs[:, end_idx].astype(\"float32\")\n    start_prefix = np.zeros((n_pairs, n_dates), dtype=\"float32\")\n    start_prefix[:, 1:] = promo_cs[:, :-1]\n    promo_fwd16_mat = (promo_fwd16_mat - start_prefix).astype(\"float32\")\n    # RAM: .copy() releases promo_fwd16 backing\n    promo_fwd16_mat = promo_fwd16_mat[:, TRAIN_COL_START:].copy()\n    print(f\"[phase 7] feature prep done \u2014 entering training loop (16 models)  \"\n          f\"RAM={psutil.Process().memory_info().rss/1e9:.2f}GB  t={time.time()-_t0:.1f}s\", flush=True)\n\n    # \u2500\u2500 3. train 16 per-lead-time models \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    # For model h (prediction day = origin + h \u2212 1):\n    #   lag-k feature = sales[origin \u2212 k] = sales_mat[:, t_idx \u2212 (h + k \u2212 1)]\n    # Constraint: t_idx \u2212 (h + max_lag \u2212 1) \u2265 0  \u2192  t_start = max(h + 71, n_dates \u2212 56)\n    models: list[lgb.LGBMRegressor | None] = []\n    max_lag_offset = max(LAGS) - 1   # 55\n\n    for h in range(1, HORIZON + 1):\n        t_start = max(h + max_lag_offset, n_dates - TRAIN_DAYS)\n        t_idx   = np.arange(t_start, n_dates)\n        n_valid = len(t_idx)\n\n        if n_valid < 2:\n            models.append(None)\n            continue\n\n        y = np.log1p(sales_mat[:, t_idx]).flatten()\n\n        lag_feats = [sales_mat[:, t_idx - (h + k - 1)].flatten() for k in LAGS]\n\n        target_dows = pivot_dates[t_idx].dayofweek.values.astype(\"float32\")\n        dow     = np.tile(target_dows, n_pairs)\n        onpromo = promo_mat[:, t_idx].flatten()\n        w       = np.repeat(pair_weights, n_valid)\n\n        # Rolling + trend features at the training origin = t_idx \u2212 (h \u2212 1)\n        origs = t_idx - (h - 1)   # (n_valid,) origin column indices\n        t_origs = origs - TRAIN_COL_START   # index into trimmed matrices\n        roll_feats  = [roll_train[name][:, t_origs].flatten() for name in _ROLL_FEAT_NAMES]\n        trend_feats = [roll_train[name][:, t_origs].flatten() for name in _TREND_FEAT_NAMES]\n        # F4: promo_share*, days_since_promo from roll_train; promo_count_fwd16 from promo_fwd16_mat\n        f4_feats = (\n            [roll_train[f\"promo_share{ps}\"][:, t_origs].flatten() for ps in _PROMO_SHARE_WINDOWS]\n            + [roll_train[\"days_since_promo\"][:, t_origs].flatten()]\n            + [promo_fwd16_mat[:, t_origs].flatten()]\n        )\n\n        # F5: Calendar features at TARGET dates (pivot_dates[t_idx])\n        f5_feats = [\n            np.tile(pivot_dom[t_idx],    n_pairs),\n            np.tile(pivot_payday[t_idx], n_pairs),\n            np.tile(pivot_days_from[t_idx], n_pairs),\n            # month/is_nat_holiday/days_to_nat_hol/is_local_holiday: bottom-decile or zero gain\n        ]\n\n        f6_feats = [hier_mats[name][:, t_origs].flatten() for name in _F6_FEAT_NAMES]\n        f7_feats = [roll_train[name][:, t_origs].flatten() for name in _F7_FEAT_NAMES]\n        f8_feats = [f8_mats[name][:, t_origs].flatten() for name in _F8_FEAT_NAMES]\n        f9_feats = [f9_lift_mat[:, t_origs].flatten()]\n\n        X = np.column_stack(lag_feats + [dow, onpromo] + roll_feats + trend_feats + f4_feats + f5_feats + f6_feats + f7_feats + f8_feats + f9_feats)\n\n        # RAM: free intermediates before fit; malloc_trim returns pages to OS\n        del lag_feats, roll_feats, trend_feats, f4_feats, f5_feats\n        del f6_feats, f7_feats, f8_feats, f9_feats, dow, onpromo\n        gc.collect(); _libc.malloc_trim(0)\n        X = X.astype(\"float32\", copy=False)   # guard: coerce any silent float64 upcast\n        if h == 1:\n            print(f\"[h=1] X dtype={X.dtype}  nbytes={X.nbytes/1e9:.2f}GB  shape={X.shape}\", flush=True)\n        # RAM logpoint (post-del, pre-fit)\n        print(f\"[h={h:2d}] RAM={psutil.Process().memory_info().rss/1e9:.2f}GB  t={time.time()-_t0:.1f}s\", flush=True)\n        clf = lgb.LGBMRegressor(**LGB_PARAMS)\n        with warnings.catch_warnings():\n            warnings.simplefilter(\"ignore\")\n            clf.fit(X, y, sample_weight=w, feature_name=FEATURE_NAMES)\n        del X, y, w; gc.collect()   # RAM: release ~4 GB before next iteration\n        models.append(clf)\n\n    # \u2500\u2500 4. inference \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n    # Lag-k at cutoff: column (n_dates \u2212 k) in the pivot.\n    lag_base = np.column_stack([sales_mat[:, n_dates - k] for k in LAGS])\n\n    # Rolling features at inference (computed from the full pivot up to cutoff)\n    roll_inf = _rolling_inf_vecs(sales_mat, promo_mat, n_dates)\n\n    # Promo lookup: EXCEPTION 3 (backtest) \u2192 EXCEPTION 1 (submission) \u2192 PROMO_NA_CODE\n    def _build_promo_lkp(df: pd.DataFrame) -> dict:\n        if df.empty:\n            return {}\n        tmp = df.copy()\n        tmp[\"promo_f\"] = _promo_float(tmp[\"onpromotion\"])\n        return tmp.set_index([\"store_nbr\", \"item_nbr\", \"date\"])[\"promo_f\"].to_dict()\n\n    train_promo_lkp = _build_promo_lkp(train_future_df)\n    test_promo_lkp  = _build_promo_lkp(test_df)\n\n    def _promo_at(store, item, date) -> float:\n        key = (store, item, date)\n        # Priority: EXCEPTION 3 (backtest) \u2192 EXCEPTION 1 (submission) \u2192 0.0.\n        # Absent pairs default to 0.0 (not PROMO_NA_CODE) so that the absence\n        # of a row in train_future_df cannot be used as a zero-sales signal.\n        v = train_promo_lkp.get(key)\n        if v is not None:\n            return float(v)\n        return float(test_promo_lkp.get(key, 0.0))\n\n    pairs_list = list(scoring_pairs.itertuples(index=False, name=None))\n\n    # F4: promo_count_fwd16 \u2014 compute once across all horizon_dates (ALLOWED EXCEPTION)\n    promo_count_fwd16_inf = np.array(\n        [sum(_promo_at(s, i, d) for d in horizon_dates) for s, i in pairs_list],\n        dtype=\"float32\",\n    )\n\n    # F4: promo_share*, days_since_promo \u2014 already in roll_inf from _rolling_inf_vecs\n    f4_inf_feats = (\n        [roll_inf[f\"promo_share{ps}\"] for ps in _PROMO_SHARE_WINDOWS]\n        + [roll_inf[\"days_since_promo\"]]\n        + [promo_count_fwd16_inf]\n    )\n\n    # F6: Hierarchy aggregates at inference \u2014 computed once (rm30 at cutoff is fixed)\n    rm30_inf = roll_inf[\"rm30\"]\n    f6_inf_feats = [\n        _hier_agg_vec(rm30_inf, grp_store),\n        _hier_agg_vec(rm30_inf, grp_item),\n        _hier_agg_vec(rm30_inf, grp_family),\n        _hier_agg_vec(rm30_inf, grp_class),\n        # cluster_agg_rm30: pruned (degenerate when 1 store/cluster)\n    ]\n\n    # F7: Sparsity at inference \u2014 computed once (all from roll_inf, fixed at cutoff)\n    f7_inf_feats = [roll_inf[name] for name in _F7_FEAT_NAMES]\n\n    # F8: External features at inference \u2014 trailing windows ending at cutoff\n    _oil_rm30_inf  = float(oil_raw[max(0, n_dates - 30) : n_dates].mean())\n    _oil_rm7_inf   = float(oil_raw[max(0, n_dates - 7)  : n_dates].mean())\n    _oil_r7_30_inf = float(np.clip(_oil_rm7_inf / (_oil_rm30_inf + 1e-9), 0.0, 5.0))\n    oil_rm30_inf_vec  = np.full(n_pairs, _oil_rm30_inf,  dtype=\"float32\")\n    oil_r7_30_inf_vec = np.full(n_pairs, _oil_r7_30_inf, dtype=\"float32\")\n    txn_rm30_inf_vec  = np.zeros(n_pairs, dtype=\"float32\")\n    txn_r7_30_inf_vec = np.zeros(n_pairs, dtype=\"float32\")\n    for _si, _s in enumerate(_ust):\n        _mask = scoring_pairs[\"store_nbr\"].values == _s\n        _row  = txn_wide.iloc[_si].values   # (n_dates,) raw txn\n        _rm30 = float(_row[max(0, n_dates - 30) : n_dates].mean())\n        _rm7  = float(_row[max(0, n_dates - 7)  : n_dates].mean())\n        txn_rm30_inf_vec[_mask]  = _rm30\n        txn_r7_30_inf_vec[_mask] = float(np.clip(_rm7 / (_rm30 + 1e-9), 0.0, 5.0))\n    f8_inf_feats = [oil_rm30_inf_vec, oil_r7_30_inf_vec, txn_rm30_inf_vec, txn_r7_30_inf_vec]\n\n    # F9: lift is constant across horizon dates \u2014 same vector used for every h\n    f9_inf_feats = [f9_lift_vec]\n\n    pred_chunks: list[pd.DataFrame] = []\n    for h, clf in enumerate(models, start=1):\n        d = horizon_dates[h - 1]\n        dow_val = float(d.dayofweek)\n\n        onpromo_inf = np.array(\n            [_promo_at(s, i, d) for s, i in pairs_list],\n            dtype=\"float32\",\n        )\n        roll_inf_feats  = [roll_inf[name] for name in _ROLL_FEAT_NAMES]\n        trend_inf_feats = [roll_inf[name] for name in _TREND_FEAT_NAMES]\n\n        # F5: Calendar features at this horizon date\n        d_date = d.date()\n        d_idx  = pd.DatetimeIndex([d])\n        _dt, _df = _nat_hol_dists(d_idx, nat_set)\n        f5_inf_feats = [\n            np.full(n_pairs, float(d.day),   \"float32\"),\n            np.full(n_pairs, _payday_prox(d_idx)[0], \"float32\"),\n            np.full(n_pairs, _df[0],  \"float32\"),\n            # month/is_nat_holiday/days_to_nat_hol/is_local_holiday: bottom-decile or zero gain\n        ]\n\n        X_inf = np.column_stack([\n            lag_base,\n            np.full(n_pairs, dow_val, dtype=\"float32\"),\n            onpromo_inf,\n        ] + roll_inf_feats + trend_inf_feats + f4_inf_feats + f5_inf_feats + f6_inf_feats + f7_inf_feats + f8_inf_feats + f9_inf_feats)\n\n        if clf is None:\n            preds = np.zeros(n_pairs, dtype=\"float32\")\n        else:\n            with warnings.catch_warnings():\n                warnings.simplefilter(\"ignore\")\n                preds = np.maximum(\n                    np.expm1(clf.predict(X_inf)), 0.0\n                ).astype(\"float32\")\n\n        pred_chunks.append(pd.DataFrame({\n            \"store_nbr\":       [s for s, _ in pairs_list],\n            \"item_nbr\":        [i for _, i in pairs_list],\n            \"date\":            d,\n            \"unit_sales_pred\": preds,\n        }))\n\n    print(f\"[done] all 16 models inferred  t={time.time()-_t0:.1f}s\", flush=True)\n    # save gain importances; set FAVORITA_IMPORTANCE_PATH=\"\" to skip\n    _fi_path = os.environ.get(\"FAVORITA_IMPORTANCE_PATH\",\n                              \"/kaggle/working/feature_importances_full.csv\")\n    if _fi_path:\n        _fi_df = feature_importances_avg(models).reset_index()\n        _fi_df.columns = [\"feature\", \"gain_avg\"]\n        _fi_df.to_csv(_fi_path, index=False)\n        print(f\"[importance] {len(_fi_df)} features \u2192 {_fi_path}\", flush=True)\n        print(_fi_df.head(20).to_string(index=False), flush=True)\n    return pd.concat(pred_chunks, ignore_index=True)\n\n\ndef feature_importances_avg(models: list) -> pd.Series:\n    \"\"\"Average gain importance across all non-None trained models.\"\"\"\n    valid = [m for m in models if m is not None]\n    if not valid:\n        return pd.Series(dtype=float)\n    fi_stack = np.vstack(\n        [m.booster_.feature_importance(importance_type=\"gain\") for m in valid]\n    )\n    return pd.Series(fi_stack.mean(axis=0), index=FEATURE_NAMES).sort_values(\n        ascending=False\n    )\n"}, {"cell_type": "markdown", "metadata": {}, "source": "## Pipeline overview\n\n**Feature engineering:** The model builds a store\u00d7item\u00d7day sales pivot from 125M+ training rows, then engineers 72 features across 9 families: lagged sales (F1); rolling mean/median/std/min/max over 5 windows 7\u2013140 d (F2); trend diffs and slopes (F3); trailing promo rate and forward promo count (F4); calendar (F5); store/item/family/class hierarchy aggregates (F6); sparsity indicators (F7); oil price and transaction volume (F8); family-level on/off-promo sales ratio (F9).\n\n**Feature selection:** 82 features (F1\u2013F9) \u2192 \u22123 structural/near-zero \u2192 79 \u2192 \u22127 bottom-decile gain ablation \u2192 **72 final**. See feature-list comments in the model cell for per-feature rationale.\n\n**Training:** 16 separate LightGBM models (one per horizon day h=1\u202616), each trained on a 56-day rolling window ending at the forecast cutoff. Three members (seeds 42, 7, 1300) are run sequentially and averaged. Blending improved NWRMSLE by ~0.0003 over a single seed in held-out backtest (2017-06-14, 2017-07-05, 2017-07-26)."}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "import sys\nsys.path.insert(0, \"/kaggle/working\")\nprint(\"sys.path[0]:\", sys.path[0])"}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "import os, time as _t\nimport pandas as pd\nfrom models.lgbm_stage3 import lgbm_stage3, LGB_PARAMS\n\n# Finale config: w2c22-pruned 72-feature\nW2C22_72F = dict(\n    objective=\"regression\", metric=\"rmse\",\n    n_estimators=200, learning_rate=0.05,\n    num_leaves=63, min_child_samples=100,\n    bagging_fraction=1.0, bagging_freq=5,\n    feature_fraction=0.8, verbose=-1, n_jobs=-1,\n)\nif SMOKE_MODE:\n    W2C22_72F[\"n_estimators\"] = 5\n    cutoff        = pd.Timestamp(\"2013-01-15\")\n    horizon_dates = pd.date_range(\"2013-01-15\", periods=16, freq=\"D\")\n    print(f\"SMOKE: n_estimators={W2C22_72F['n_estimators']}  cutoff={cutoff.date()}\")\nelse:\n    cutoff        = pd.Timestamp(\"2017-08-16\")\n    horizon_dates = pd.date_range(\"2017-08-16\", periods=16, freq=\"D\")\n    print(f\"cutoff: {cutoff.date()}  \"\n          f\"horizon: {horizon_dates[0].date()} \u2192 {horizon_dates[-1].date()}\")\n\nt_train_start = time.time()\n\n# Seeds were fixed before this submission; not tuned on the public leaderboard.\n# Blending 3 seeds improved NWRMSLE by ~0.0003 vs single-seed in held-out backtest\n# (origins 2017-06-14, 2017-07-05, 2017-07-26).\n\n# Member A \u2014 seed=42, saves feature importances\nLGB_PARAMS.update({**W2C22_72F, \"seed\": 42})\nos.environ[\"FAVORITA_IMPORTANCE_PATH\"] = \"/kaggle/working/feature_importances_full.csv\"\nt0 = _t.time()\npreds_A = lgbm_stage3(cutoff, horizon_dates, dev=SMOKE_MODE)\nprint(f\"[blend] member A (seed=42):   {_t.time()-t0:.0f}s  shape={preds_A.shape}\")\nos.environ[\"FAVORITA_IMPORTANCE_PATH\"] = \"\"\n\n# Member B \u2014 seed=7\nLGB_PARAMS.update({\"seed\": 7})\nt0 = _t.time()\npreds_B = lgbm_stage3(cutoff, horizon_dates, dev=SMOKE_MODE)\nprint(f\"[blend] member B (seed=7):    {_t.time()-t0:.0f}s\")\n\n# Member C \u2014 seed=1300\nLGB_PARAMS.update({\"seed\": 1300})\nt0 = _t.time()\npreds_C = lgbm_stage3(cutoff, horizon_dates, dev=SMOKE_MODE)\nprint(f\"[blend] member C (seed=1300): {_t.time()-t0:.0f}s\")\n\n# Blend \u2014 simple mean\npreds = preds_A.copy()\npreds[\"unit_sales_pred\"] = (\n    preds_A[\"unit_sales_pred\"] +\n    preds_B[\"unit_sales_pred\"] +\n    preds_C[\"unit_sales_pred\"]\n) / 3\nprint(f\"\\n[blend] {len(preds):,} rows  total={time.time()-t_train_start:.0f}s\")\nprint(f\"Training + blend wall-clock: {time.time() - t_train_start:.1f}s\")\nprint(f\"Predictions shape:   {preds.shape}\")\n"}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "from dataguard import load\n\n# allow_future=True: test.csv includes the onpromo column for the forecast window.\n# The competition provides it \u2014 promotion status is known at prediction time by design.\n# This is intentional (competition Exception 4), not a data-leak workaround.\ntest_manifest = load(\"test\", dev=False, allow_future=True)[\n    [\"id\", \"store_nbr\", \"item_nbr\", \"date\"]\n]\nsub = test_manifest.merge(\n    preds[[\"store_nbr\", \"item_nbr\", \"date\", \"unit_sales_pred\"]],\n    on=[\"store_nbr\", \"item_nbr\", \"date\"],\n    how=\"left\",\n)\nsub[\"unit_sales\"] = sub[\"unit_sales_pred\"].fillna(0.0).clip(lower=0.0)\nassert len(sub) == 3_370_464, f\"Row count mismatch: {len(sub)}\"\nsub[[\"id\", \"unit_sales\"]].to_csv(\"submission.csv\", index=False)\nprint(f\"Submission rows: {len(sub):,}\")\nprint(\"submission.csv written.\")"}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "print(f\"Total notebook wall-clock: {time.time() - t_notebook_start:.1f}s\")\nprint(f\"n_pairs: {len(preds[['store_nbr','item_nbr']].drop_duplicates()):,}\")\n\nif SMOKE_MODE:\n    print()\n    print(\"=\" * 60)\n    print(\"  SMOKE MODE COMPLETE \u2014 all cells executed successfully\")\n    print(\"=\" * 60)"}]}