{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b37ff653",
   "metadata": {},
   "source": [
    "# Workspace price-tier prediction: final notebook\n",
    "\n",
    "**Task.** Each row is a bookable workspace listing. Predict its `price_tier` (1\u20134). Metric: **macro F1** (the average F1 of the 4 tiers, so every tier counts equally).\n",
    "\n",
    "**The key finding.** Kaggle split train and test **by operator**: only ~1% of test rows belong to an operator that appears in train.\n",
    "Our first model validated with random K-fold. That puts the same operator in both training and validation, so the model got\n",
    "rewarded for memorising operators. It scored **0.758 in CV but 0.70 on the leaderboard**.\n",
    "Validating with folds **grouped by operator** reproduces the leaderboard (\u22480.70), so every decision here is made on grouped CV.\n",
    "\n",
    "**Pipeline**\n",
    "1. Load data, rebuild an operator id, make operator-grouped folds\n",
    "2. EDA: target, metros, size, terms, and the operator split\n",
    "3. Preprocessing & features: encodings, ratios, geo, label-free operator profile, out-of-fold target features\n",
    "4. Feature-block experiments (measured earlier, summarised)\n",
    "5. Models: Logistic Regression (scaled pipeline baseline), XGBoost, XGBoost regression, CatBoost (+ LightGBM optional)\n",
    "6. Level-2 stacker: smooth predictions across each operator's listings\n",
    "7. Scenarios: blends \u00d7 decoders, all cross-fitted on grouped CV; the best one is submitted\n",
    "8. Submission\n",
    "9. Error analysis\n",
    "10. What else could be tried\n",
    "\n",
    "**How to run.** Colab: Runtime \u2192 Change runtime type \u2192 **T4 GPU**, then Run all (the first cell fetches the data). Kaggle: attach the competition data,\n",
    "set Accelerator \u2192 GPU T4, then Run all. Every finished fold is checkpointed in `ckpt/`, and a rerun skips work that's already done."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bca9e82a",
   "metadata": {},
   "source": [
    "## 0. Setup (Colab only; skipped on Kaggle)\n",
    "Runtime \u2192 Change runtime type \u2192 **T4 GPU** first. This cell installs CatBoost and downloads the competition data from the course Drive link."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21258cbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys, subprocess, zipfile\n",
    "from pathlib import Path as _P\n",
    "if \"google.colab\" in sys.modules:\n",
    "    subprocess.run([sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"catboost\", \"gdown\"], check=True)\n",
    "    if not any(_P(\"/content\").rglob(\"train.csv\")):\n",
    "        import gdown\n",
    "        gdown.download(id=\"1v6Ce6HLFvByXzDdcHQD1NMgA-XnJf5kf\", output=\"/content/data_download\", quiet=False)\n",
    "        z = _P(\"/content/data_download\")\n",
    "        if zipfile.is_zipfile(z):\n",
    "            zipfile.ZipFile(z).extractall(\"/content/data\")\n",
    "            for inner in list(_P(\"/content/data\").rglob(\"*.zip\")):  # zip inside a zip\n",
    "                zipfile.ZipFile(inner).extractall(inner.parent)\n",
    "    print(\"train.csv found at:\", [str(p) for p in _P(\"/content\").rglob(\"train.csv\")])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07795218",
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "import os\n",
    "import shutil\n",
    "import subprocess\n",
    "import time\n",
    "import warnings\n",
    "from pathlib import Path\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import lightgbm as lgb\n",
    "import xgboost as xgb\n",
    "from catboost import CatBoostClassifier\n",
    "from scipy import sparse\n",
    "from scipy.optimize import minimize\n",
    "from sklearn.cluster import MiniBatchKMeans\n",
    "from sklearn.metrics import log_loss\n",
    "from sklearn.model_selection import StratifiedGroupKFold\n",
    "from sklearn.neighbors import NearestNeighbors\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "# ======================= SETTINGS =======================\n",
    "QUICK_TEST = bool(int(os.environ.get(\"PT_QUICK\", \"0\")))  # True = 10% of operators + few rounds (pipeline check only)\n",
    "USE_LIGHTGBM = False                                        # LightGBM runs on CPU only (its GPU mode is slower on Kaggle): ~8 min/fold\n",
    "RUN_ABLATION = False                                        # test H1/H2 (feature blocks) before the big training\n",
    "FOLDS, SEED, K = 5, 42, 4\n",
    "MAX_ROUNDS = 300 if QUICK_TEST else 6000\n",
    "LOG_EVERY = 100 if QUICK_TEST else 250                     # print training progress every N rounds\n",
    "# =========================================================\n",
    "\n",
    "TAG = f\"final{'q' if QUICK_TEST else ''}_f{FOLDS}\"          # checkpoint names (quick-test checkpoints never mix with full ones)\n",
    "WORK = Path(\"/kaggle/working\") if Path(\"/kaggle/working\").exists() else Path(\".\")\n",
    "CKPT = WORK / \"ckpt\"\n",
    "CKPT.mkdir(exist_ok=True)\n",
    "SEARCH = [CKPT] + (sorted(p for p in Path(\"/kaggle/input\").glob(\"**/ckpt\") if p.is_dir()) if Path(\"/kaggle/input\").exists() else [])\n",
    "\n",
    "\n",
    "def gpu_available():\n",
    "    if shutil.which(\"nvidia-smi\") is None:\n",
    "        return False\n",
    "    try:\n",
    "        return subprocess.run([\"nvidia-smi\"], capture_output=True, timeout=10).returncode == 0\n",
    "    except Exception:\n",
    "        return False\n",
    "\n",
    "\n",
    "GPU = gpu_available()\n",
    "if GPU:  # all-GPU models; LightGBM only if asked for\n",
    "    MODELS = [\"logreg\", \"xgb\", \"xgb_reg\", \"cat\"] + ([\"lgb_multi\", \"lgb_reg\"] if USE_LIGHTGBM else [])\n",
    "else:\n",
    "    MODELS = [\"logreg\", \"lgb_multi\", \"lgb_reg\", \"xgb\", \"xgb_reg\"] + ([\"cat\"] if os.environ.get(\"PT_FORCE_CAT\") else [])\n",
    "    print(\"No GPU found: running on CPU and skipping CatBoost. Turn on GPU T4 for a much faster run.\")\n",
    "N_THREADS = os.cpu_count()\n",
    "print(f\"GPU: {GPU} | CPU threads: {N_THREADS} | models: {MODELS} | quick test: {QUICK_TEST}\")\n",
    "print(\"checkpoint folders searched:\", [str(p) for p in SEARCH])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c62bf211",
   "metadata": {},
   "source": [
    "## Helpers: checkpoints, progress bar, fast macro F1\n",
    "Checkpoints are plain `.npz` files (NumPy arrays only, no pickle). They are written to a temp file and then renamed,\n",
    "so a crash during a save never leaves a broken checkpoint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5236de3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ck_load(name):\n",
    "    for d in SEARCH:\n",
    "        p = d / f\"{name}.npz\"\n",
    "        if p.exists():\n",
    "            with np.load(p, allow_pickle=False) as z:\n",
    "                return {k: z[k] for k in z.files}\n",
    "    return None\n",
    "\n",
    "\n",
    "def ck_save(name, **arrays):\n",
    "    tmp = CKPT / f\"{name}.tmp.npz\"\n",
    "    np.savez(tmp, **arrays)\n",
    "    tmp.replace(CKPT / f\"{name}.npz\")\n",
    "\n",
    "\n",
    "class Progress:\n",
    "    \"\"\"Overall progress bar over all big steps. ETA = remaining steps x the measured time of the same kind of step\n",
    "    (e.g. each remaining XGBoost fold is estimated from finished XGBoost folds, not from the slower LightGBM ones).\"\"\"\n",
    "\n",
    "    def __init__(self, plan):\n",
    "        self.plan, self.total, self.done, self.t0, self.kind_times = plan, len(plan), 0, time.time(), {}\n",
    "\n",
    "    def bar(self, width=30):\n",
    "        full = int(width * self.done / self.total)\n",
    "        return \"\u2588\" * full + \"\u2591\" * (width - full) + f\" {self.done}/{self.total} ({100 * self.done / self.total:.0f}%)\"\n",
    "\n",
    "    def start(self, label):\n",
    "        self.label, self.ts = label, time.time()\n",
    "        print(f\"\\n\u25b6 {label}   [{self.bar()}]\", flush=True)\n",
    "\n",
    "    def finish(self, from_checkpoint=False):\n",
    "        self.done += 1\n",
    "        took = time.time() - self.ts\n",
    "        if not from_checkpoint:\n",
    "            self.kind_times.setdefault(self.label.split(\":\")[0], []).append(took)\n",
    "        left = [s.split(\":\")[0] for s in self.plan[self.done:]]\n",
    "        timed = sum(np.mean(self.kind_times[k]) for k in left if k in self.kind_times)\n",
    "        untimed = sum(k not in self.kind_times for k in left)\n",
    "        eta = \"-\" if not left else f\"{timed / 60:.0f} min\" + (f\" + {untimed} steps not timed yet\" if untimed else \"\")\n",
    "        how = \"loaded from checkpoint\" if from_checkpoint else f\"took {took / 60:.1f} min\"\n",
    "        print(f\"\u2713 {self.label}: {how} | [{self.bar()}] | elapsed {(time.time() - self.t0) / 60:.1f} min | ETA ~{eta}\", flush=True)\n",
    "\n",
    "\n",
    "def conf(yt, yp, w=None):\n",
    "    return np.bincount(yt * K + yp, weights=w, minlength=K * K).reshape(K, K)\n",
    "\n",
    "\n",
    "def f1_from_conf(cm):\n",
    "    tp = np.diag(cm)\n",
    "    return float(np.mean(2 * tp / np.maximum(cm.sum(0) + cm.sum(1), 1e-12)))\n",
    "\n",
    "\n",
    "def fast_f1(yt, yp, w=None):\n",
    "    \"\"\"macro F1 (optionally sample-weighted), ~50\u00d7 faster than sklearn; used inside the tuning loops\"\"\"\n",
    "    return f1_from_conf(conf(yt, yp, w))\n",
    "\n",
    "\n",
    "PLAN = ([\"load data + grouped folds\", \"build features\", \"ablation of feature blocks\"]\n",
    "        + [f\"{m}: fold {f + 1}/{FOLDS}\" for m in MODELS for f in range(FOLDS)]\n",
    "        + [\"level-2 operator stacker\", \"scenario table\", \"write submissions\"])\n",
    "prog = Progress(PLAN)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39f220c9",
   "metadata": {},
   "source": [
    "## 1. Data, operator id, grouped folds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e41b4cae",
   "metadata": {},
   "outputs": [],
   "source": [
    "prog.start(\"load data + grouped folds\")\n",
    "\n",
    "\n",
    "DATA_DIR = \"\"  # optional: paste the folder that contains train.csv, e.g. \"/kaggle/input/workstation-rental-price-tier-prediction\"\n",
    "\n",
    "\n",
    "def find_data():\n",
    "    if DATA_DIR:\n",
    "        return Path(DATA_DIR)\n",
    "    if \"PT_DATA\" in os.environ:\n",
    "        return Path(os.environ[\"PT_DATA\"])\n",
    "    for root in [\"/kaggle/input\", \"/content\", \".\"]:\n",
    "        for dirpath, _, files in os.walk(root, followlinks=True):  # followlinks: Kaggle mounts inputs as symlinks\n",
    "            if \"train.csv\" in files and \"test.csv\" in files:\n",
    "                return Path(dirpath)\n",
    "    found = [str(q) for q in Path(\"/kaggle/input\").glob(\"*\")] if Path(\"/kaggle/input\").exists() else []\n",
    "    raise FileNotFoundError(\n",
    "        \"train.csv / test.csv not found. The competition data is not attached to this notebook. \"\n",
    "        \"Fix: right panel -> 'Add Input' -> 'Competitions' -> search 'workstation-rental-price-tier-prediction' -> Add. \"\n",
    "        f\"Currently in /kaggle/input: {found}\")\n",
    "\n",
    "\n",
    "DATA = find_data()\n",
    "train = pd.read_csv(DATA / \"train.csv\")\n",
    "test = pd.read_csv(DATA / \"test.csv\")\n",
    "\n",
    "# operator_ref alone isn't unique (different operators share names); these columns together identify one operator\n",
    "HOST = [\"operator_ref\", \"operator_platform_age_years\", \"operator_platform_age_months\", \"operator_active_years\",\n",
    "        \"operator_active_months\", \"operator_class\", \"operator_portfolio_size\"]\n",
    "host = pd.concat([train[HOST], test[HOST]], ignore_index=True).groupby(HOST, dropna=False, sort=False).ngroup().to_numpy()\n",
    "if QUICK_TEST:  # keep 10% of operators (whole operators, so the split logic stays the same)\n",
    "    keep = (host * 2654435761 % 1000) < 100\n",
    "    n0 = len(train)\n",
    "    train, test = train[keep[:n0]].reset_index(drop=True), test[keep[n0:]].reset_index(drop=True)\n",
    "    host = host[keep]\n",
    "\n",
    "n_tr = len(train)\n",
    "y = train[\"price_tier\"].to_numpy() - 1\n",
    "folds = list(StratifiedGroupKFold(FOLDS, shuffle=True, random_state=SEED).split(np.zeros(n_tr), y, host[:n_tr]))\n",
    "fold_id = np.zeros(n_tr, int)\n",
    "for i, (_, v) in enumerate(folds):\n",
    "    fold_id[v] = i\n",
    "print(train.shape, test.shape, \"| test rows whose operator is in train:\", round(np.isin(host[n_tr:], host[:n_tr]).mean(), 4))\n",
    "prog.finish()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e88b897",
   "metadata": {},
   "source": [
    "## 2. Exploratory data analysis\n",
    "Goal: understand the target, the strongest drivers, data quality, and **how train and test differ** before building anything."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8db53002",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "pd.set_option(\"display.width\", 160)\n",
    "print(\"rows:\", len(train), \"train /\", len(test), \"test | columns:\", train.shape[1])\n",
    "print(train.dtypes.astype(str).value_counts().to_dict())\n",
    "miss = train.isna().mean().sort_values(ascending=False)\n",
    "display(miss[miss > 0].round(3).to_frame(\"share missing\").head(15))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8afcaf3",
   "metadata": {},
   "source": [
    "**Missing values** are concentrated in review dates/scores (listings never reviewed) and a few size columns. They're informative\n",
    "(no reviews = new or unpopular listing), so tree models keep them as \"missing\" rather than filling in a guess."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc1b1699",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 3, figsize=(18, 4))\n",
    "train.price_tier.value_counts().sort_index().plot.bar(ax=ax[0], title=\"price_tier (target)\", rot=0)\n",
    "pd.crosstab(train.access_type, train.price_tier, normalize=\"index\").plot.barh(stacked=True, ax=ax[1], title=\"tier mix by access type\")\n",
    "sns.boxplot(data=train.assign(seats=train.max_seats.clip(upper=16)), x=\"price_tier\", y=\"seats\", ax=ax[2], showfliers=False)\n",
    "ax[2].set_title(\"max_seats by tier\"); plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aed003d3",
   "metadata": {},
   "source": [
    "- The 4 tiers are roughly balanced (tier 2 is the largest), so **macro F1** treats them all equally and no class can be ignored.\n",
    "- **Access type** (whole space vs private vs shared zone) shifts the tier mix strongly.\n",
    "- **Size** (`max_seats`) climbs steadily with the tier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95d0c9fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(16, 5))\n",
    "pd.crosstab(train.metro, train.price_tier, normalize=\"index\").plot.barh(stacked=True, ax=ax[0], title=\"tier mix by metro\", legend=False)\n",
    "shift = pd.DataFrame({\"train\": train.metro.value_counts(normalize=True), \"test\": test.metro.value_counts(normalize=True)})\n",
    "shift.sort_values(\"train\").plot.barh(ax=ax[1], title=\"metro share: train vs test\")\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c73e7361",
   "metadata": {},
   "source": [
    "- **Metro matters enormously**: some metros are mostly tier 1, others mostly tier 4 (tiers are price bands, and price levels differ by market).\n",
    "- The **metro mix differs between train and test** (e.g. one metro is ~4% of train but ~10% of test), a first hint that the split isn't random."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c5bf618",
   "metadata": {},
   "outputs": [],
   "source": [
    "num = [\"max_seats\", \"zones\", \"workstations\", \"utility_rooms\", \"min_term_days\", \"booked_load\", \"review_count\",\n",
    "       \"score_overall\", \"score_value\", \"operator_portfolio_size\", \"open_slots_d\"]\n",
    "corr = train[num + [\"price_tier\"]].corr(method=\"spearman\")[\"price_tier\"].drop(\"price_tier\").sort_values()\n",
    "corr.plot.barh(figsize=(7, 4), title=\"Spearman correlation with price_tier\"); plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd952e92",
   "metadata": {},
   "source": [
    "Size columns (seats, zones, workstations, utility rooms) correlate most with the tier. Long **minimum terms** (monthly rentals) go with cheaper tiers.\n",
    "No single column is decisive. The model has to combine size, location, terms and the operator type."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdc8f812",
   "metadata": {},
   "source": [
    "### The operator split: the most important finding\n",
    "Listings are run by **operators**, and an operator's listings tend to share a tier. How did Kaggle split train/test with respect to operators?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8522f3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import StratifiedKFold\n",
    "g_ = pd.Series(y).groupby(host[:n_tr])\n",
    "big = g_.size() >= 5\n",
    "print(f\"operators with \u22655 listings: share of their listings in the operator's most common tier = \"\n",
    "      f\"{g_.agg(lambda s: s.value_counts(normalize=True).iloc[0])[big].mean():.2f}\")\n",
    "rnd = list(StratifiedKFold(FOLDS, shuffle=True, random_state=SEED).split(np.zeros(n_tr), y))\n",
    "seen_rnd = np.mean([np.isin(host[v], host[f]).mean() for f, v in rnd])\n",
    "seen_grp = np.mean([np.isin(host[v], host[f]).mean() for f, v in folds])\n",
    "seen_test = np.isin(host[n_tr:], host[:n_tr]).mean()\n",
    "display(pd.DataFrame({\"share of rows whose operator was seen in training\": [seen_rnd, seen_grp, seen_test]},\n",
    "                     index=[\"random K-fold validation\", \"operator-grouped validation\", \"real test set\"]).round(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1789e6d",
   "metadata": {},
   "source": [
    "**This explains everything about our first model:**\n",
    "- Operators are very consistent (about 3 in 4 of their listings share one tier).\n",
    "- With **random folds**, most validation rows have their operator in training, so the model can \"recognise\" the operator. That's an easy exam.\n",
    "- In the **real test set**, almost no operator was seen in training. That's the hard exam.\n",
    "\n",
    "So random-fold CV (0.758) was measuring the wrong thing. **Grouped folds** (an operator is entirely in train or entirely in validation) match the test set,\n",
    "and gave \u22480.70, the same as the leaderboard. From here on, every score is grouped CV."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef98aab6",
   "metadata": {},
   "source": [
    "## 3. Preprocessing & features (checkpointed)\n",
    "**Preprocessing decisions**\n",
    "- **Missing values**: tree models (XGBoost, CatBoost, LightGBM) learn which side of a split missing values go to, so we keep NaN.\n",
    "  The Logistic Regression baseline uses a pipeline with `SimpleImputer(median)`.\n",
    "- **Scaling**: trees only use the *order* of values, so they need no scaling. Logistic Regression does (`StandardScaler` in its pipeline).\n",
    "- **Encoding**: low-cardinality text columns (metro, access type, space class, \u2026) \u2192 integer codes; high-cardinality ones (operator, zone, borough) \u2192 frequency counts\n",
    "  plus **out-of-fold target encodings** (the tier mix of the zone / comparable listings), computed on grouped folds so no row ever sees its own label.\n",
    "- Everything is built once for train+test with the same functions, so both get identical transformations. Target statistics come from training folds only.\n",
    "\n",
    "Blocks (the ablation in section 3 decides which to keep):\n",
    "- **v2 set**: base columns, facility one-hots (\u2265 0.5% of listings), geo transforms, demand/space ratios, local density, zone context,\n",
    "  label-free operator profile, and out-of-fold target features on grouped folds (zone / comparables / geo neighbours / kNN comparables).\n",
    "- **H1 `rel_*`**: listing minus the average of its comparables (same zone + access type + seats), for bookings, open days, reviews, terms, facilities, scores. No labels used.\n",
    "- **H2 `facte_*`**: out-of-fold facility effects. For each facility, the average tier residual after removing the metro average; per listing, the mean / sum / max / min over its facilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9dd0bf2",
   "metadata": {},
   "outputs": [],
   "source": [
    "prog.start(\"build features\")\n",
    "feat = ck_load(f\"{TAG}_features\")\n",
    "if feat is None:\n",
    "    raw = pd.concat([train.drop(columns=\"price_tier\"), test], ignore_index=True)\n",
    "    rawx = raw.assign(seats_b=raw.max_seats.clip(upper=12))\n",
    "    df = raw.copy()\n",
    "    for c in [\"operator_certified\", \"operator_has_avatar\", \"operator_id_verified\", \"bookable\"]:\n",
    "        df[c] = df[c].map({\"t\": 1, \"f\": 0})\n",
    "    for c in [\"snapshot_start_date\", \"snapshot_end_date\", \"first_review_date\", \"last_review_date\"]:\n",
    "        df[c] = (pd.to_datetime(df[c]) - pd.Timestamp(\"2000-01-01\")).dt.days\n",
    "    for c in [\"operator_ref\", \"zone\", \"borough\"]:\n",
    "        df[c] = raw[c].map(raw[c].value_counts())\n",
    "    nfac = raw[\"facilities\"].str.count(\"FCL-\").fillna(0)\n",
    "    df[\"n_facilities\"] = nfac\n",
    "    df[\"occupancy\"] = df.booked_load / 365\n",
    "    df[\"review_speed_12m\"] = df.review_count_12m / (df.review_count + 1)\n",
    "    df[\"review_trend\"] = (df.review_count_12m + 1) / (df.review_count_prev_yr + 1)\n",
    "    df[\"days_since_last_review\"] = df.snapshot_end_date - df.last_review_date\n",
    "    df[\"listing_age_days\"] = df.snapshot_end_date - df.first_review_date\n",
    "    df[\"open_ratio_a_d\"] = df.open_slots_a / (df.open_slots_d + 1)\n",
    "    df[\"open_ratio_e_d\"] = df.open_slots_e / (df.open_slots_d + 1)\n",
    "    df[\"seats_per_ws\"] = df.max_seats / df.workstations.replace(0, np.nan)\n",
    "    df[\"seats_per_zone\"] = df.max_seats / df.zones.replace(0, np.nan)\n",
    "    df[\"seats_per_util\"] = df.max_seats / df.utility_rooms.replace(0, np.nan)\n",
    "    df[\"term_span\"] = df.max_term_days - df.min_term_days\n",
    "\n",
    "    codes = raw[\"facilities\"].fillna(\"[]\").map(json.loads).explode().dropna()\n",
    "    fac = pd.crosstab(codes.index, codes).reindex(raw.index, fill_value=0).clip(upper=1).astype(np.uint8)\n",
    "    fac = fac.loc[:, fac.mean() >= 0.005]\n",
    "    fac.columns = [\"f_\" + c[4:] for c in fac.columns]\n",
    "    code_id, code_names = pd.factorize(codes.values)\n",
    "    F_all = sparse.csr_matrix((np.ones(len(codes), np.float32), (codes.index.to_numpy(), code_id)), shape=(len(raw), len(code_names)))\n",
    "    F_all.data[:] = 1.0\n",
    "    F_all = F_all[:, np.where(np.asarray(F_all.sum(0)).ravel() >= 30)[0]].tocsr()\n",
    "\n",
    "    lat, lon = raw[\"geo_lat\"], raw[\"geo_lon\"]\n",
    "    geo = pd.DataFrame({\"lat_p_lon\": lat + lon, \"lat_m_lon\": lat - lon})\n",
    "    for a in (15, 30, 60):\n",
    "        r = np.deg2rad(a)\n",
    "        geo[f\"rot{a}_x\"] = lat * np.cos(r) - lon * np.sin(r)\n",
    "        geo[f\"rot{a}_y\"] = lat * np.sin(r) + lon * np.cos(r)\n",
    "    for k in (\"metro\", \"zone\"):\n",
    "        center = raw.groupby(k)[[\"geo_lat\", \"geo_lon\"]].transform(\"median\")\n",
    "        geo[f\"dist_{k}_center\"] = np.hypot(lat - center[\"geo_lat\"], lon - center[\"geo_lon\"])\n",
    "    km = MiniBatchKMeans(100, random_state=0, n_init=3, batch_size=4096).fit(raw[[\"geo_lat\", \"geo_lon\"]])\n",
    "    geo[\"geo_cluster\"] = km.labels_\n",
    "    geo[\"geo_cluster_dist\"] = km.transform(raw[[\"geo_lat\", \"geo_lon\"]]).min(1)\n",
    "    dist_all, _ = NearestNeighbors(n_neighbors=51).fit(raw[[\"geo_lat\", \"geo_lon\"]]).kneighbors(raw[[\"geo_lat\", \"geo_lon\"]])\n",
    "    geo[\"density_d10\"], geo[\"density_d50\"] = dist_all[:, 10], dist_all[:, 50]\n",
    "    geo[\"zone_mean_seats\"] = df.max_seats.groupby(raw.zone).transform(\"mean\")\n",
    "    geo[\"zone_mean_nfac\"] = nfac.groupby(raw.zone).transform(\"mean\")\n",
    "    geo[\"seats_vs_zone\"] = df.max_seats - geo.zone_mean_seats\n",
    "    geo[\"nfac_vs_zone\"] = nfac - geo.zone_mean_nfac\n",
    "\n",
    "    # label-free operator profile: what kind of operator is this, and how does this listing compare with its portfolio?\n",
    "    d = pd.DataFrame({\"seats\": raw.max_seats, \"ws\": raw.workstations, \"zones\": raw.zones, \"util\": raw.utility_rooms,\n",
    "                      \"nfac\": nfac, \"score\": raw.score_overall, \"rev\": raw.review_count,\n",
    "                      \"minterm\": raw.min_term_days, \"open_d\": raw.open_slots_d, \"load\": raw.booked_load})\n",
    "    for c in [\"access_type\", \"space_class\"]:\n",
    "        for v in raw[c].value_counts().index[:4]:\n",
    "            d[f\"{c}_{v}\"] = (raw[c] == v).astype(float)\n",
    "    g = d.groupby(host)\n",
    "    prof = pd.DataFrame(index=d.index)\n",
    "    prof[\"h_rows\"] = g.seats.transform(\"size\")\n",
    "    prof[\"h_zones_n\"] = raw.groupby(host).zone.transform(\"nunique\")\n",
    "    prof[\"h_metros_n\"] = raw.groupby(host).metro.transform(\"nunique\")\n",
    "    for c in d.columns:\n",
    "        m = g[c].transform(\"mean\")\n",
    "        prof[f\"h_mean_{c}\"] = m\n",
    "        if c in (\"seats\", \"ws\", \"zones\", \"util\", \"nfac\", \"score\", \"minterm\"):\n",
    "            prof[f\"h_std_{c}\"] = g[c].transform(\"std\")\n",
    "            prof[f\"h_diff_{c}\"] = d[c] - m\n",
    "    prof[\"h_rank_seats\"] = g.seats.rank(pct=True)\n",
    "    prof[\"h_rank_nfac\"] = g.nfac.rank(pct=True)\n",
    "    bld = raw.groupby([\"geo_lat\", \"geo_lon\"]).ngroup().to_numpy()\n",
    "    gb = d.groupby(bld)\n",
    "    prof[\"b_rows\"] = gb.seats.transform(\"size\")\n",
    "    prof[\"b_hosts\"] = pd.Series(host).groupby(bld).transform(\"nunique\")\n",
    "    prof[\"b_mean_seats\"] = gb.seats.transform(\"mean\")\n",
    "    prof[\"b_mean_nfac\"] = gb.nfac.transform(\"mean\")\n",
    "\n",
    "    # H1: demand and quality relative to comparable listings (same zone + access type + seats); no labels\n",
    "    comp = rawx.groupby([\"zone\", \"access_type\", \"seats_b\"], dropna=False, sort=False).ngroup().to_numpy()\n",
    "    value_gap = raw.score_value - raw[[\"score_accuracy\", \"score_condition\", \"score_onboarding\", \"score_support\", \"score_access\"]].mean(axis=1)\n",
    "    rel = pd.DataFrame({\"rel_value_gap\": value_gap, \"rel_comp_size\": pd.Series(comp).map(pd.Series(comp).value_counts()).to_numpy()})\n",
    "    for c, s in {\"load\": raw.booked_load, \"open_a\": raw.open_slots_a, \"open_d\": raw.open_slots_d, \"rpm\": raw.reviews_per_month,\n",
    "                 \"rev12\": raw.review_count_12m, \"minterm\": raw.min_term_days, \"nfac\": nfac, \"score\": raw.score_overall,\n",
    "                 \"value\": value_gap, \"ws\": raw.workstations}.items():\n",
    "        rel[f\"rel_{c}_vs_comps\"] = s - s.groupby(comp).transform(\"mean\")\n",
    "    metro_comp = rawx.groupby([\"metro\", \"access_type\", \"seats_b\"], dropna=False, sort=False).ngroup().to_numpy()\n",
    "    rel[\"rel_load_vs_metro_comps\"] = raw.booked_load - raw.booked_load.groupby(metro_comp).transform(\"mean\")\n",
    "\n",
    "    df = df.drop(columns=[\"facilities\", \"space_id\"])\n",
    "    CAT_COLS = list(df.select_dtypes(exclude=\"number\").columns)\n",
    "    for c in CAT_COLS:\n",
    "        df[c] = df[c].astype(\"category\").cat.codes  # -1 = missing (LightGBM treats negative categories as missing)\n",
    "    base = pd.concat([df, fac, geo, prof, rel], axis=1)\n",
    "\n",
    "    # ---- target features, out-of-fold on grouped folds ----\n",
    "    TE_GROUPS = {\n",
    "        \"zone\": [\"zone\"], \"boro\": [\"borough\"], \"sc\": [\"space_class\"], \"metro_sc\": [\"metro\", \"space_class\"],\n",
    "        \"zone_acc\": [\"zone\", \"access_type\"], \"zone_seats\": [\"zone\", \"max_seats\"], \"zone_sc\": [\"zone\", \"space_class\"],\n",
    "        \"geo\": [\"geo_lat\", \"geo_lon\"], \"metro_acc_seats\": [\"metro\", \"access_type\", \"max_seats\"],\n",
    "        \"zone_acc_seats\": [\"zone\", \"access_type\", \"seats_b\"],\n",
    "        \"zone_acc_seats_ur\": [\"zone\", \"access_type\", \"seats_b\", \"utility_rooms\"],\n",
    "        \"boro_acc_seats\": [\"borough\", \"access_type\", \"seats_b\"], \"metro_sc_seats\": [\"metro\", \"space_class\", \"seats_b\"],\n",
    "        \"bld_acc\": [\"geo_lat\", \"geo_lon\", \"access_type\"], \"metro_acc_rooms\": [\"metro\", \"access_type\", \"zones\", \"utility_rooms\"],\n",
    "    }\n",
    "    RES_GROUPS = [\"zone\", \"geo\", \"zone_sc\", \"sc\"]\n",
    "    gid = {k: rawx[c].groupby(c, sort=False, dropna=False).ngroup().to_numpy() for k, c in TE_GROUPS.items()}\n",
    "    base_ids = raw[[\"access_type\", \"max_seats\"]].groupby([\"access_type\", \"max_seats\"], dropna=False).ngroup().to_numpy()\n",
    "    metro_id = pd.factorize(raw.metro)[0]\n",
    "    coords = raw[[\"geo_lat\", \"geo_lon\"]].to_numpy()\n",
    "    prior = np.bincount(y, minlength=K) / len(y)\n",
    "    tiers = np.arange(1, K + 1)\n",
    "    space = np.column_stack([(lat - lat.mean()) / lat.std(), (lon - lon.mean()) / lon.std(), np.log1p(raw.max_seats) * 0.3,\n",
    "                             raw.zones.fillna(1) * 0.05, raw.utility_rooms.fillna(1) * 0.05, nfac * 0.002]).astype(np.float32)\n",
    "    acc = raw.access_type.astype(str).to_numpy()\n",
    "\n",
    "    def target_features(fit, apply):\n",
    "        cols = {}\n",
    "        for name, g_ in gid.items():  # smoothed tier distribution of each group\n",
    "            cnt = np.zeros((g_.max() + 1, K)); np.add.at(cnt, (g_[fit], y[fit]), 1)\n",
    "            tot = cnt.sum(1)\n",
    "            p = ((cnt + 20 * prior) / (tot[:, None] + 20))[g_[apply]]\n",
    "            for j in range(K):\n",
    "                cols[f\"te_{name}_p{j + 1}\"] = p[:, j]\n",
    "            cols[f\"te_{name}_mean\"] = p @ tiers\n",
    "            cols[f\"te_{name}_n\"] = np.log1p(tot[g_[apply]])\n",
    "        yt = y[fit] + 1.0  # residuals vs. an (access type, seats) baseline\n",
    "        nb = base_ids.max() + 1\n",
    "        bm = (np.bincount(base_ids[fit], yt, nb) + 20 * yt.mean()) / (np.bincount(base_ids[fit], minlength=nb) + 20)\n",
    "        res = yt - bm[base_ids[fit]]\n",
    "        for name in RES_GROUPS:\n",
    "            g_ = gid[name]; n = g_.max() + 1\n",
    "            c_ = np.bincount(g_[fit], minlength=n)[g_[apply]] + 10\n",
    "            mres = np.bincount(g_[fit], res, n)[g_[apply]] / c_\n",
    "            cols[f\"res_{name}_mean\"] = mres\n",
    "            cols[f\"res_{name}_pred\"] = bm[base_ids[apply]] + mres\n",
    "        pts, inv = np.unique(coords[fit], axis=0, return_inverse=True)  # tiers of nearest locations\n",
    "        cnt = np.zeros((len(pts), K)); np.add.at(cnt, (inv.ravel(), y[fit]), 1)\n",
    "        dist, idx = NearestNeighbors(n_neighbors=50).fit(pts).kneighbors(coords[apply])\n",
    "        for k in (1, 5, 20, 50):\n",
    "            c_ = cnt[idx[:, :k]].sum(1)\n",
    "            p = (c_ + 5 * prior) / (c_.sum(1)[:, None] + 5)\n",
    "            for j in range(K):\n",
    "                cols[f\"geo{k}_p{j + 1}\"] = p[:, j]\n",
    "            cols[f\"geo{k}_mean\"] = p @ tiers\n",
    "            cols[f\"geo{k}_dist\"] = dist[:, k - 1]\n",
    "        for k in (10, 50):  # comparable listings: nearest in (location, size) with the same access type\n",
    "            for s in (\"mean\", \"p1\", \"p4\", \"dist\"):\n",
    "                cols[f\"knn{k}_{s}\"] = np.full(len(apply), np.nan)\n",
    "        for a in np.unique(acc):\n",
    "            fa, aa = fit[acc[fit] == a], np.where(acc[apply] == a)[0]\n",
    "            if len(fa) < 50 or len(aa) == 0:\n",
    "                continue\n",
    "            dist, idx = NearestNeighbors(n_neighbors=50).fit(space[fa]).kneighbors(space[apply[aa]])\n",
    "            lab = y[fa][idx]\n",
    "            for k in (10, 50):\n",
    "                oh = np.stack([(lab[:, :k] == j).mean(1) for j in range(K)], 1)\n",
    "                cols[f\"knn{k}_mean\"][aa] = oh @ tiers\n",
    "                cols[f\"knn{k}_p1\"][aa], cols[f\"knn{k}_p4\"][aa] = oh[:, 0], oh[:, 3]\n",
    "                cols[f\"knn{k}_dist\"][aa] = dist[:, k - 1]\n",
    "        # H2: facility effects after removing the metro average (shrunk towards 0 for rare facilities)\n",
    "        nm = metro_id.max() + 1\n",
    "        mm = np.bincount(metro_id[fit], yt, nm) / np.maximum(np.bincount(metro_id[fit], minlength=nm), 1)\n",
    "        r_ = yt - mm[metro_id[fit]]\n",
    "        Ff = F_all[fit]\n",
    "        dev = (Ff.T @ r_) / (np.asarray(Ff.sum(0)).ravel() + 50)\n",
    "        Fa = F_all[apply]\n",
    "        nf = np.asarray(Fa.sum(1)).ravel()\n",
    "        has = nf > 0\n",
    "        cols[\"facte_sum\"] = Fa @ dev\n",
    "        cols[\"facte_mean\"] = np.where(has, cols[\"facte_sum\"] / np.maximum(nf, 1), np.nan)\n",
    "        cols[\"facte_max\"] = np.where(has, Fa.multiply(dev + 10).max(1).toarray().ravel() - 10, np.nan)\n",
    "        cols[\"facte_min\"] = np.where(has, 10 - Fa.multiply(10 - dev).max(1).toarray().ravel(), np.nan)\n",
    "        cols[\"facte_n_premium\"] = Fa @ (dev > 0.05).astype(np.float32)\n",
    "        cols[\"facte_n_budget\"] = Fa @ (dev < -0.05).astype(np.float32)\n",
    "        return pd.DataFrame(cols)\n",
    "\n",
    "    T_tr = pd.concat([target_features(f, v).set_index(v) for f, v in folds]).sort_index()\n",
    "    T_te = target_features(np.arange(n_tr), np.arange(n_tr, len(raw)))\n",
    "    X_df = pd.concat([base.iloc[:n_tr].reset_index(drop=True), T_tr], axis=1)\n",
    "    Xt_df = pd.concat([base.iloc[n_tr:].reset_index(drop=True), T_te], axis=1)\n",
    "    feat = dict(X=X_df.to_numpy(np.float32), Xt=Xt_df.to_numpy(np.float32), cols=np.array(X_df.columns, dtype=str),\n",
    "                cat_cols=np.array(CAT_COLS, dtype=str), bld=bld)\n",
    "    ck_save(f\"{TAG}_features\", **feat)\n",
    "    del raw, rawx, df, base, X_df, Xt_df, T_tr, T_te, F_all\n",
    "    prog.finish()\n",
    "else:\n",
    "    prog.finish(from_checkpoint=True)\n",
    "\n",
    "X_full, Xt_full = feat[\"X\"], feat[\"Xt\"]\n",
    "ALL_COLS = list(feat[\"cols\"])\n",
    "bld = feat[\"bld\"]\n",
    "MT_tr = np.maximum(X_full[:, ALL_COLS.index(\"metro\")], 0).astype(int)\n",
    "MT_te = np.maximum(Xt_full[:, ALL_COLS.index(\"metro\")], 0).astype(int)\n",
    "print(\"all features:\", X_full.shape, \"| test:\", Xt_full.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19381844",
   "metadata": {},
   "source": [
    "## 4. Feature-block experiments\n",
    "Two extra feature blocks were tested with a fast model on 50% of operators (grouped CV, full data, measured on Kaggle):\n",
    "\n",
    "| feature set | features | grouped CV macro F1 |\n",
    "|---|---|---|\n",
    "| v2 (base + geo + operator profile + target features) | 392 | **0.7028** |\n",
    "| v2 + H1 demand relative to comparables (`rel_*`) | 405 | 0.7016 |\n",
    "| v2 + H2 facility price effects (`facte_*`) | 398 | 0.7029 |\n",
    "| v2 + both | 411 | 0.7026 |\n",
    "\n",
    "Neither block beats noise (\u00b10.002), so the simpler **v2** set is used (`RUN_ABLATION = False` skips re-running this).\n",
    "Lesson: a hypothesis that is plausible and even visible in the data can still add nothing once the model already has related features.\n",
    "\n",
    "Fast settings on 50% of operators and 3 folds (XGBoost on GPU, otherwise LightGBM).\n",
    "Each result is checkpointed. Differences under ~0.002 are within noise, and then the smaller set is kept."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "368b99f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "prog.start(\"ablation of feature blocks\")\n",
    "FEATURE_SETS = {\n",
    "    \"v2\": [c for c in ALL_COLS if not c.startswith((\"rel_\", \"facte_\"))],\n",
    "    \"v2 + H1 relative\": [c for c in ALL_COLS if not c.startswith(\"facte_\")],\n",
    "    \"v2 + H2 facility\": [c for c in ALL_COLS if not c.startswith(\"rel_\")],\n",
    "    \"v2 + H1 + H2\": ALL_COLS,\n",
    "}\n",
    "abl_rows = []\n",
    "if RUN_ABLATION:\n",
    "    sub = (host[:n_tr] * 2654435761 % 1000) < 500\n",
    "    for name, cols in FEATURE_SETS.items():\n",
    "        idx = [ALL_COLS.index(c) for c in cols]\n",
    "        scores = []\n",
    "        for f in range(3):\n",
    "            ckn = f\"{TAG}_abl_{name.replace(' ', '_').replace('+', 'p')}_{f}\"\n",
    "            ck = ck_load(ckn)\n",
    "            if ck is None:\n",
    "                tr = np.where(sub & (fold_id != f))[0]\n",
    "                va = np.where(sub & (fold_id == f))[0]\n",
    "                if GPU:\n",
    "                    pa = dict(objective=\"multi:softprob\", num_class=K, eta=0.1, max_depth=8, min_child_weight=5, subsample=0.8,\n",
    "                              colsample_bytree=0.5, max_bin=128, tree_method=\"hist\", device=\"cuda\", eval_metric=\"mlogloss\", seed=SEED)\n",
    "                    dva = xgb.DMatrix(X_full[va][:, idx], label=y[va])\n",
    "                    b = xgb.train(pa, xgb.DMatrix(X_full[tr][:, idx], label=y[tr]), 3000, evals=[(dva, \"v\")],\n",
    "                                  early_stopping_rounds=50, verbose_eval=False)\n",
    "                    p = b.predict(dva, iteration_range=(0, b.best_iteration + 1))\n",
    "                else:\n",
    "                    pa = dict(objective=\"multiclass\", num_class=K, learning_rate=0.1, num_leaves=127, min_child_samples=40,\n",
    "                              feature_fraction=0.5, bagging_fraction=0.8, bagging_freq=1, max_bin=63, verbose=-1, seed=SEED, num_threads=N_THREADS)\n",
    "                    m = lgb.train(pa, lgb.Dataset(X_full[tr][:, idx], y[tr]), 3000, valid_sets=[lgb.Dataset(X_full[va][:, idx], y[va])],\n",
    "                                  callbacks=[lgb.early_stopping(50, verbose=False)])\n",
    "                    p = m.predict(X_full[va][:, idx], num_iteration=m.best_iteration)\n",
    "                ck = {\"f1\": np.array([fast_f1(y[va], p.argmax(1))])}\n",
    "                ck_save(ckn, **ck)\n",
    "            scores.append(float(ck[\"f1\"][0]))\n",
    "        abl_rows.append({\"feature set\": name, \"n features\": len(cols), \"grouped CV macro F1\": round(np.mean(scores), 4)})\n",
    "        print(abl_rows[-1], flush=True)\n",
    "    abl = pd.DataFrame(abl_rows)\n",
    "    best = abl[\"grouped CV macro F1\"].max()\n",
    "    # simplest set within 0.002 of the best (fewer features = faster + less overfitting)\n",
    "    FINAL_SET = abl[abl[\"grouped CV macro F1\"] >= best - 0.002].sort_values(\"n features\").iloc[0][\"feature set\"]\n",
    "    display(abl)\n",
    "else:\n",
    "    FINAL_SET = \"v2\"  # winner of the measured experiments above\n",
    "print(\"chosen feature set:\", FINAL_SET)\n",
    "SEL = [ALL_COLS.index(c) for c in FEATURE_SETS[FINAL_SET]]\n",
    "COLS = FEATURE_SETS[FINAL_SET]\n",
    "X, Xt = X_full[:, SEL], Xt_full[:, SEL]\n",
    "CAT_IDX = [COLS.index(c) for c in feat[\"cat_cols\"] if c in COLS]\n",
    "SET_TAG = FINAL_SET.replace(\" \", \"\").replace(\"+\", \"p\")\n",
    "prog.finish()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ac484b6",
   "metadata": {},
   "source": [
    "## 5. Models (each fold checkpointed)\n",
    "| model | output | why it's here |\n",
    "|---|---|---|\n",
    "| `logreg` Logistic Regression (impute \u2192 scale \u2192 LR pipeline) | 4 probabilities | simple, fast, interpretable **baseline**; shows what non-linear models add |\n",
    "| `lgb_multi` LightGBM multiclass | 4 probabilities | strong default for tabular data |\n",
    "| `lgb_reg` LightGBM regression on tier 1\u20134 | one score | uses the **order** of the tiers (H4) |\n",
    "| `xgb` XGBoost multiclass | 4 probabilities | a different boosting implementation, so its mistakes differ |\n",
    "| `xgb_reg` XGBoost regression on tier 1\u20134 | one score | a second ordered model, on GPU |\n",
    "| `cat` CatBoost multiclass | 4 probabilities | symmetric trees, strong regularisation, a third opinion |\n",
    "| `lgb_ord` (optional) | 4 probabilities | v1's ordinal model: 3 binary models P(tier > k) |\n",
    "\n",
    "Every model uses early stopping on its validation fold; test predictions are averaged over the 5 fold models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5133fd0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def pick_lgb_device():\n",
    "    \"\"\"LightGBM's GPU mode (OpenCL) isn't always available or faster: time 30 rounds on 60k rows on each device, keep the faster.\"\"\"\n",
    "    if not GPU or not any(m.startswith(\"lgb\") for m in MODELS):\n",
    "        return \"cpu\"\n",
    "    rows_ = np.random.RandomState(0).choice(n_tr, min(60000, n_tr), replace=False)\n",
    "    times = {}\n",
    "    for dev in (\"cpu\", \"gpu\"):\n",
    "        try:\n",
    "            pa = dict(objective=\"multiclass\", num_class=K, learning_rate=0.05, num_leaves=127, max_bin=63, feature_fraction=0.5,\n",
    "                      verbose=-1, num_threads=N_THREADS, device_type=dev)\n",
    "            t0 = time.time()\n",
    "            lgb.train(pa, lgb.Dataset(X[rows_], y[rows_], categorical_feature=CAT_IDX), 30)\n",
    "            times[dev] = time.time() - t0\n",
    "        except Exception as e:\n",
    "            print(f\"LightGBM on {dev} not available: {str(e)[:150]}\")\n",
    "    print(\"LightGBM speed test (30 rounds, 60k rows):\", {k: f\"{v:.1f}s\" for k, v in times.items()})\n",
    "    return min(times, key=times.get) if times else \"cpu\"\n",
    "\n",
    "\n",
    "LGB_DEVICE = pick_lgb_device()\n",
    "print(\"LightGBM device:\", LGB_DEVICE, \"| XGBoost device:\", \"cuda\" if GPU else \"cpu\", \"| CatBoost:\", \"GPU\" if GPU else \"CPU\")\n",
    "LGB_BASE = dict(learning_rate=0.05, num_leaves=127, min_child_samples=40, feature_fraction=0.5, bagging_fraction=0.8,\n",
    "                bagging_freq=1, lambda_l2=1.0, max_bin=63, verbose=-1, seed=SEED, num_threads=N_THREADS, device_type=LGB_DEVICE)\n",
    "XGB_PARAMS = dict(objective=\"multi:softprob\", num_class=K, eta=0.05, max_depth=8, min_child_weight=5, subsample=0.8,\n",
    "                  colsample_bytree=0.5, max_bin=128, tree_method=\"hist\", eval_metric=\"mlogloss\", seed=SEED,\n",
    "                  device=\"cuda\" if GPU else \"cpu\", nthread=N_THREADS)\n",
    "CAT_PARAMS = dict(loss_function=\"MultiClass\", iterations=MAX_ROUNDS, learning_rate=0.08, depth=8, border_count=128,\n",
    "                  random_seed=SEED, task_type=\"GPU\" if GPU else \"CPU\", thread_count=N_THREADS, od_type=\"Iter\", od_wait=100)\n",
    "\n",
    "\n",
    "def train_lgb(params, Xtr, ytr, Xva, yva, Xte):\n",
    "    dtr = lgb.Dataset(Xtr, ytr, categorical_feature=CAT_IDX)\n",
    "    dva = lgb.Dataset(Xva, yva, categorical_feature=CAT_IDX, reference=dtr)\n",
    "    m = lgb.train(params, dtr, MAX_ROUNDS, valid_sets=[dva],\n",
    "                  callbacks=[lgb.early_stopping(100, verbose=False), lgb.log_evaluation(LOG_EVERY)])\n",
    "    print(f\"   best iteration: {m.best_iteration}\")\n",
    "    return m.predict(Xva, num_iteration=m.best_iteration), m.predict(Xte, num_iteration=m.best_iteration)\n",
    "\n",
    "\n",
    "def cumulative_to_classes(p_above):\n",
    "    p_above = np.minimum.accumulate(p_above, axis=1)\n",
    "    edges = np.column_stack([np.ones(len(p_above)), p_above, np.zeros(len(p_above))])\n",
    "    return np.clip(edges[:, :-1] - edges[:, 1:], 1e-6, 1)\n",
    "\n",
    "\n",
    "def fit_lgb_multi(Xtr, ytr, Xva, yva, Xte):\n",
    "    return train_lgb({**LGB_BASE, \"objective\": \"multiclass\", \"num_class\": K}, Xtr, ytr, Xva, yva, Xte)\n",
    "\n",
    "\n",
    "def fit_lgb_reg(Xtr, ytr, Xva, yva, Xte):\n",
    "    pv, pt = train_lgb({**LGB_BASE, \"objective\": \"regression\"}, Xtr, ytr + 1.0, Xva, yva + 1.0, Xte)\n",
    "    return pv[:, None], pt[:, None]\n",
    "\n",
    "\n",
    "def fit_lgb_ord(Xtr, ytr, Xva, yva, Xte):\n",
    "    av, at = [], []\n",
    "    for k in range(K - 1):\n",
    "        v, t = train_lgb({**LGB_BASE, \"objective\": \"binary\"}, Xtr, (ytr > k).astype(int), Xva, (yva > k).astype(int), Xte)\n",
    "        av.append(v); at.append(t)\n",
    "    return cumulative_to_classes(np.column_stack(av)), cumulative_to_classes(np.column_stack(at))\n",
    "\n",
    "\n",
    "from sklearn.impute import SimpleImputer\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "\n",
    "def fit_logreg(Xtr, ytr, Xva, yva, Xte):\n",
    "    \"\"\"baseline: impute -> scale -> multinomial logistic regression (a reusable sklearn pipeline)\"\"\"\n",
    "    clean = lambda A: np.where(np.isfinite(A), A, np.nan)\n",
    "    m = make_pipeline(SimpleImputer(strategy=\"median\"), StandardScaler(), LogisticRegression(max_iter=300, C=0.5))\n",
    "    m.fit(clean(Xtr), ytr)\n",
    "    return m.predict_proba(clean(Xva)), m.predict_proba(clean(Xte))\n",
    "\n",
    "\n",
    "LAST_IMP = {}\n",
    "\n",
    "\n",
    "def fit_xgb(Xtr, ytr, Xva, yva, Xte):\n",
    "    dtr, dva = xgb.DMatrix(Xtr, label=ytr), xgb.DMatrix(Xva, label=yva)\n",
    "    b = xgb.train(XGB_PARAMS, dtr, MAX_ROUNDS, evals=[(dva, \"valid\")], early_stopping_rounds=100, verbose_eval=LOG_EVERY)\n",
    "    LAST_IMP.clear(); LAST_IMP.update(b.get_score(importance_type=\"gain\"))\n",
    "    it = (0, b.best_iteration + 1)\n",
    "    print(f\"   best iteration: {b.best_iteration}\")\n",
    "    return b.predict(dva, iteration_range=it), b.predict(xgb.DMatrix(Xte), iteration_range=it)\n",
    "\n",
    "\n",
    "XGB_REG = {**{k: v for k, v in XGB_PARAMS.items() if k not in (\"objective\", \"num_class\", \"eval_metric\")},\n",
    "           \"objective\": \"reg:squarederror\", \"eval_metric\": \"rmse\"}\n",
    "\n",
    "\n",
    "def fit_xgb_reg(Xtr, ytr, Xva, yva, Xte):\n",
    "    dtr, dva = xgb.DMatrix(Xtr, label=ytr + 1.0), xgb.DMatrix(Xva, label=yva + 1.0)\n",
    "    b = xgb.train(XGB_REG, dtr, MAX_ROUNDS, evals=[(dva, \"valid\")], early_stopping_rounds=100, verbose_eval=LOG_EVERY)\n",
    "    it = (0, b.best_iteration + 1)\n",
    "    print(f\"   best iteration: {b.best_iteration}\")\n",
    "    return b.predict(dva, iteration_range=it)[:, None], b.predict(xgb.DMatrix(Xte), iteration_range=it)[:, None]\n",
    "\n",
    "\n",
    "def fit_cat(Xtr, ytr, Xva, yva, Xte):\n",
    "    m = CatBoostClassifier(**CAT_PARAMS)\n",
    "    m.fit(Xtr, ytr, eval_set=(Xva, yva), verbose=LOG_EVERY, use_best_model=True)\n",
    "    return m.predict_proba(Xva), m.predict_proba(Xte)\n",
    "\n",
    "\n",
    "TRAINERS = {\"logreg\": fit_logreg, \"lgb_multi\": fit_lgb_multi, \"lgb_reg\": fit_lgb_reg, \"lgb_ord\": fit_lgb_ord, \"xgb\": fit_xgb, \"xgb_reg\": fit_xgb_reg,\n",
    "            \"cat\": fit_cat}\n",
    "OOF, TEST = {}, {}\n",
    "for key in MODELS:\n",
    "    for f, (tr, va) in enumerate(folds):\n",
    "        name = f\"{TAG}_{SET_TAG}_{key}_fold{f}\"\n",
    "        prog.start(f\"{key}: fold {f + 1}/{FOLDS}\")\n",
    "        ck = ck_load(name)\n",
    "        if ck is None:\n",
    "            pv, pt = TRAINERS[key](X[tr], y[tr], X[va], y[va], Xt)\n",
    "            ck = {\"oof\": np.asarray(pv, np.float32), \"test\": np.asarray(pt, np.float32)}\n",
    "            if key == \"xgb\":  # keep feature importance (gain) for section 5c\n",
    "                ck[\"imp\"] = np.array([LAST_IMP.get(f\"f{i}\", 0.0) for i in range(X.shape[1])], np.float32)\n",
    "            ck_save(name, **ck)\n",
    "            prog.finish()\n",
    "        else:\n",
    "            prog.finish(from_checkpoint=True)\n",
    "        if key not in OOF:\n",
    "            OOF[key] = np.zeros((n_tr, ck[\"oof\"].shape[1]), np.float32)\n",
    "            TEST[key] = np.zeros((len(Xt), ck[\"test\"].shape[1]), np.float32)\n",
    "        OOF[key][va] = ck[\"oof\"]\n",
    "        if \"imp\" in ck:\n",
    "            XGB_IMP = XGB_IMP + ck[\"imp\"] if \"XGB_IMP\" in globals() else ck[\"imp\"].copy()\n",
    "        TEST[key] += ck[\"test\"] / FOLDS\n",
    "        if ck[\"oof\"].shape[1] == K:\n",
    "            print(f\"   fold macro F1: {fast_f1(y[va], ck['oof'].argmax(1)):.4f}\")\n",
    "for key in OOF:\n",
    "    if OOF[key].shape[1] == K:\n",
    "        print(f\"{key:10s} grouped CV macro F1 (plain argmax): {fast_f1(y, OOF[key].argmax(1)):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb7e90e6",
   "metadata": {},
   "source": [
    "### 5b. Model comparison\n",
    "Every model is scored on the same operator-grouped folds. Probability models use the most likely tier; regression models round to the nearest tier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f8f374c",
   "metadata": {},
   "outputs": [],
   "source": [
    "cmp_rows = []\n",
    "for key in OOF:\n",
    "    pred_ = OOF[key].argmax(1) if OOF[key].shape[1] == K else np.clip(np.rint(OOF[key][:, 0]) - 1, 0, K - 1).astype(int)\n",
    "    cmp_rows.append({\"model\": key, \"grouped CV macro F1\": round(fast_f1(y, pred_), 4)})\n",
    "display(pd.DataFrame(cmp_rows).sort_values(\"grouped CV macro F1\", ascending=False).reset_index(drop=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fac185e",
   "metadata": {},
   "source": [
    "| model | advantages | disadvantages |\n",
    "|---|---|---|\n",
    "| Logistic Regression | fast, stable, coefficients are interpretable, little overfitting | only linear boundaries; needs scaling/imputation; can't use interactions like \"zone \u00d7 size\" unless we build them |\n",
    "| XGBoost | strong on tabular data, handles missing values, captures interactions, fast on GPU | many hyperparameters; can overfit without early stopping |\n",
    "| XGBoost regression | uses the **order** of tiers (predicting 2.6 is \"between 2 and 3\") | needs cut-points to turn a number into a tier |\n",
    "| CatBoost | symmetric trees regularise well; robust with default settings | slower on CPU; less flexible trees |\n",
    "| LightGBM (optional) | very fast on CPU with many cores, leaf-wise growth | its GPU mode was slower than CPU here; on Kaggle's 4 cores it is the bottleneck |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92b7fe19",
   "metadata": {},
   "source": [
    "### 5c. Which features matter? (XGBoost gain, summed over folds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49558cc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "if \"XGB_IMP\" in globals():\n",
    "    imp = pd.Series(XGB_IMP, index=COLS).sort_values(ascending=False)\n",
    "    fig, ax = plt.subplots(figsize=(8, 7))\n",
    "    imp.head(25)[::-1].plot.barh(ax=ax, title=\"Top 25 features by XGBoost gain\")\n",
    "    plt.tight_layout(); plt.show()\n",
    "    groups = {\"target encodings (te_/res_)\": (\"te_\", \"res_\"), \"geo neighbours (geo*/knn*)\": (\"geo\", \"knn\"),\n",
    "              \"operator profile (h_/b_)\": (\"h_\", \"b_\"), \"facilities (f_)\": (\"f_\",)}\n",
    "    share = {g: imp[[c for c in imp.index if c.startswith(p)]].sum() / imp.sum() for g, p in groups.items()}\n",
    "    share[\"listing's own columns\"] = 1 - sum(share.values())\n",
    "    display(pd.Series(share, name=\"share of total gain\").round(3).sort_values(ascending=False).to_frame())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e5803b8",
   "metadata": {},
   "source": [
    "## 6. Level-2 stacker: operator smoothing\n",
    "For each listing, we average the level-1 predictions of the **same operator's other listings** (leave-one-out),\n",
    "and of other listings at the **same coordinates**. A small LightGBM learns how much to trust these, per metro and access type.\n",
    "With grouped folds every operator sits entirely inside one validation fold, so this is measured honestly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d153e2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "prog.start(\"level-2 operator stacker\")\n",
    "PROB_MODELS = [k for k in OOF if OOF[k].shape[1] == K]\n",
    "REG_MODELS = [k for k in OOF if OOF[k].shape[1] == 1]\n",
    "\n",
    "\n",
    "def loo_mean(P, key):\n",
    "    _, inv = np.unique(key, return_inverse=True)\n",
    "    S = np.zeros((inv.max() + 1, P.shape[1])); np.add.at(S, inv, P)\n",
    "    other = np.bincount(inv)[inv] - 1\n",
    "    M = (S[inv] - P) / np.maximum(other, 1)[:, None]\n",
    "    M[other == 0] = np.nan\n",
    "    return M, other\n",
    "\n",
    "\n",
    "def stack_inputs(store, host_ids, bld_ids, Xsrc):\n",
    "    parts = [store[k] for k in PROB_MODELS + REG_MODELS]\n",
    "    avg = np.mean([store[k] for k in PROB_MODELS], axis=0)\n",
    "    summary = np.column_stack([avg] + [store[k] for k in REG_MODELS])\n",
    "    H, hn = loo_mean(summary, host_ids)\n",
    "    B, bn = loo_mean(summary, bld_ids)\n",
    "    ctx = Xsrc[:, [ALL_COLS.index(c) for c in (\"metro\", \"access_type\", \"max_seats\")]]\n",
    "    return np.column_stack(parts + [H, hn, np.log1p(hn), B, bn, ctx]).astype(np.float32)\n",
    "\n",
    "\n",
    "Z = stack_inputs(OOF, host[:n_tr], bld[:n_tr] * FOLDS + fold_id, X_full)  # buildings: only rows predicted together\n",
    "Zt = stack_inputs(TEST, host[n_tr:], bld[n_tr:], Xt_full)\n",
    "stack_name = f\"{TAG}_{SET_TAG}_stack_{'-'.join(MODELS)}\"\n",
    "ck = ck_load(stack_name)\n",
    "if ck is None:\n",
    "    P2 = dict(objective=\"multiclass\", num_class=K, learning_rate=0.05, num_leaves=31, min_child_samples=200,\n",
    "              feature_fraction=0.9, bagging_fraction=0.8, bagging_freq=1, verbose=-1, seed=SEED, num_threads=N_THREADS)\n",
    "    oof2, test2 = np.zeros((n_tr, K)), np.zeros((len(Zt), K))\n",
    "    for f, v in folds:\n",
    "        m = lgb.train(P2, lgb.Dataset(Z[f], y[f]), 3000, valid_sets=[lgb.Dataset(Z[v], y[v])],\n",
    "                      callbacks=[lgb.early_stopping(100, verbose=False)])\n",
    "        oof2[v] = m.predict(Z[v], num_iteration=m.best_iteration)\n",
    "        test2 += m.predict(Zt, num_iteration=m.best_iteration) / FOLDS\n",
    "    ck = {\"oof\": oof2, \"test\": test2}\n",
    "    ck_save(stack_name, **ck)\n",
    "    prog.finish()\n",
    "else:\n",
    "    prog.finish(from_checkpoint=True)\n",
    "OOF[\"stack\"], TEST[\"stack\"] = ck[\"oof\"], ck[\"test\"]\n",
    "print(f\"stack      grouped CV macro F1 (plain argmax): {fast_f1(y, OOF['stack'].argmax(1)):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07f5970d",
   "metadata": {},
   "source": [
    "## 7. Scenarios: every combination, honestly scored\n",
    "**Outputs** (H3, H5): single models, a plain average, a weight-optimised average (weights minimise log-loss), probabilities mixed with the regression score, and the stacker.\n",
    "\n",
    "**Decoders** (how an output becomes a tier):\n",
    "- **A. class weights**: multiply the 4 probabilities by tuned weights and take the largest (v1's method),\n",
    "- **B. cut-points** (H4): expected tier \u03a3 p\u00b7tier (or the regression score), cut at 3 tuned thresholds,\n",
    "- **C. per-metro cut-points** (H7): one set of thresholds per metro (tried on the 4 best B scenarios).\n",
    "\n",
    "**Scoring** (H6): every decoder is tuned with sample weights that make the train metro mix look like the test metro mix.\n",
    "It is then **cross-fitted**: tuned on 4 folds' out-of-fold predictions and scored on the 5th, so no scenario wins by tuning on its own answers.\n",
    "The main column, \"test-weighted\", estimates the score on the test set's metro mix. \"plain\" is the unweighted version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82c576e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "prog.start(\"scenario table\")\n",
    "TIERS = np.arange(1, K + 1)\n",
    "share_tr = np.bincount(MT_tr, minlength=MT_tr.max() + 1) / n_tr\n",
    "share_te = np.bincount(MT_te, minlength=MT_tr.max() + 1) / len(MT_te)\n",
    "W_tr = (share_te / np.maximum(share_tr, 1e-12))[MT_tr]\n",
    "print(\"metro weights (test share / train share):\", np.round(share_te / np.maximum(share_tr, 1e-12), 2))\n",
    "\n",
    "\n",
    "def fit_cw(P, yy, w, mt):\n",
    "    W = np.ones(K)\n",
    "    for _ in range(3):\n",
    "        for k in range(1, K):\n",
    "            grid = W[k] * np.exp(np.linspace(-0.4, 0.4, 41))\n",
    "            sc = [fast_f1(yy, (P * np.r_[W[:k], g_, W[k + 1:]]).argmax(1), w) for g_ in grid]\n",
    "            W[k] = grid[int(np.argmax(sc))]\n",
    "    return W\n",
    "\n",
    "\n",
    "def apply_cw(P, W, mt):\n",
    "    return (P * W).argmax(1)\n",
    "\n",
    "\n",
    "def fit_cuts(s, yy, w, mt=None):\n",
    "    cuts = np.quantile(s, np.cumsum(np.bincount(yy, minlength=K) / len(yy))[:-1])\n",
    "    for _ in range(3):\n",
    "        for i in range(K - 1):\n",
    "            lo = cuts[i - 1] + 1e-6 if i else s.min()\n",
    "            hi = cuts[i + 1] - 1e-6 if i < K - 2 else s.max()\n",
    "            grid = np.linspace(max(lo, cuts[i] - 0.4), min(hi, cuts[i] + 0.4), 61)\n",
    "            sc = [fast_f1(yy, np.digitize(s, np.r_[cuts[:i], g_, cuts[i + 1:]]), w) for g_ in grid]\n",
    "            cuts[i] = grid[int(np.argmax(sc))]\n",
    "    return cuts\n",
    "\n",
    "\n",
    "def apply_cuts(s, cuts, mt):\n",
    "    return np.digitize(s, cuts)\n",
    "\n",
    "\n",
    "def fit_metro_cuts(s, yy, w, mt):\n",
    "    g0 = fit_cuts(s, yy, w)\n",
    "    C = np.tile(g0, (MT_tr.max() + 1, 1))\n",
    "    pred = np.digitize(s, g0)\n",
    "    cm = conf(yy, pred, w)\n",
    "    for _ in range(2):\n",
    "        for m_ in np.unique(mt):\n",
    "            ix = mt == m_\n",
    "            for i in range(K - 1):\n",
    "                rest = cm - conf(yy[ix], pred[ix], w[ix])\n",
    "                lo = C[m_, i - 1] + 1e-6 if i else s.min()\n",
    "                hi = C[m_, i + 1] - 1e-6 if i < K - 2 else s.max()\n",
    "                grid = np.linspace(max(lo, C[m_, i] - 0.5), min(hi, C[m_, i] + 0.5), 41)\n",
    "                best_g, best_s = C[m_, i], -1\n",
    "                for g_ in grid:\n",
    "                    p_ = np.digitize(s[ix], np.r_[C[m_, :i], g_, C[m_, i + 1:]])\n",
    "                    sc = f1_from_conf(rest + conf(yy[ix], p_, w[ix]))\n",
    "                    if sc > best_s:\n",
    "                        best_g, best_s = g_, sc\n",
    "                C[m_, i] = best_g\n",
    "                pred[ix] = np.digitize(s[ix], C[m_])\n",
    "                cm = rest + conf(yy[ix], pred[ix], w[ix])\n",
    "    return C\n",
    "\n",
    "\n",
    "def apply_metro_cuts(s, C, mt):\n",
    "    out = np.empty(len(s), int)\n",
    "    for m_ in np.unique(mt):\n",
    "        ix = mt == m_\n",
    "        out[ix] = np.digitize(s[ix], C[m_])\n",
    "    return out\n",
    "\n",
    "\n",
    "def crossfit(o, fit, apply):\n",
    "    pred = np.empty(n_tr, int)\n",
    "    for f in range(FOLDS):\n",
    "        m_ = fold_id == f\n",
    "        prm = fit(o[~m_], y[~m_], W_tr[~m_], MT_tr[~m_])\n",
    "        pred[m_] = apply(o[m_], prm, MT_tr[m_])\n",
    "    return fast_f1(y, pred, W_tr), fast_f1(y, pred)\n",
    "\n",
    "\n",
    "def opt_weights(models):\n",
    "    Ps = [OOF[m] for m in models]\n",
    "    def loss(z):\n",
    "        w = np.exp(z) / np.exp(z).sum()\n",
    "        return log_loss(y, np.clip(sum(wi * p for wi, p in zip(w, Ps)), 1e-7, 1))\n",
    "    z = minimize(loss, np.zeros(len(Ps)), method=\"Nelder-Mead\", options={\"maxiter\": 300}).x\n",
    "    return np.exp(z) / np.exp(z).sum()\n",
    "\n",
    "\n",
    "candidates = {}  # name -> (oof output, test output, kind)\n",
    "for m_ in PROB_MODELS + [\"stack\"]:\n",
    "    candidates[m_] = (OOF[m_], TEST[m_], \"prob\")\n",
    "for m_ in REG_MODELS:\n",
    "    candidates[m_] = (OOF[m_][:, 0], TEST[m_][:, 0], \"score\")\n",
    "if len(PROB_MODELS) > 1:\n",
    "    candidates[\"avg(\" + \"+\".join(PROB_MODELS) + \")\"] = (np.mean([OOF[m_] for m_ in PROB_MODELS], 0),\n",
    "                                                       np.mean([TEST[m_] for m_ in PROB_MODELS], 0), \"prob\")\n",
    "    wb = opt_weights(PROB_MODELS)\n",
    "    print(\"optimised blend weights:\", dict(zip(PROB_MODELS, wb.round(3))))\n",
    "    candidates[\"opt_blend\"] = (sum(wi * OOF[m_] for wi, m_ in zip(wb, PROB_MODELS)),\n",
    "                               sum(wi * TEST[m_] for wi, m_ in zip(wb, PROB_MODELS)), \"prob\")\n",
    "if REG_MODELS:\n",
    "    reg_o = np.mean([OOF[m_][:, 0] for m_ in REG_MODELS], 0)\n",
    "    reg_t = np.mean([TEST[m_][:, 0] for m_ in REG_MODELS], 0)\n",
    "    if len(REG_MODELS) > 1:\n",
    "        candidates[\"avg(\" + \"+\".join(REG_MODELS) + \")\"] = (reg_o, reg_t, \"score\")\n",
    "    reg_opts = {m_: (OOF[m_][:, 0], TEST[m_][:, 0]) for m_ in REG_MODELS}\n",
    "    if len(REG_MODELS) > 1:\n",
    "        reg_opts[\"avg_reg\"] = (reg_o, reg_t)\n",
    "    for src in [c for c in (\"opt_blend\", \"stack\") if c in candidates]:\n",
    "        for rname, (ro, rt) in reg_opts.items():\n",
    "            for a in (0.3, 0.5, 0.7):\n",
    "                candidates[f\"{src} + {rname} (a={a})\"] = (a * (candidates[src][0] @ TIERS) + (1 - a) * ro,\n",
    "                                                          a * (candidates[src][1] @ TIERS) + (1 - a) * rt, \"score\")\n",
    "\n",
    "rows = []\n",
    "for name, (o, t, kind) in candidates.items():\n",
    "    s_o, s_t = (o @ TIERS, t @ TIERS) if kind == \"prob\" else (o, t)\n",
    "    if kind == \"prob\":\n",
    "        rows.append(dict(scenario=name, decoder=\"A class weights\", out=(o, t), fit=fit_cw, apply=apply_cw,\n",
    "                         score=crossfit(o, fit_cw, apply_cw)))\n",
    "    rows.append(dict(scenario=name, decoder=\"B cut-points\", out=(s_o, s_t), fit=fit_cuts, apply=apply_cuts,\n",
    "                     score=crossfit(s_o, fit_cuts, apply_cuts)))\n",
    "    print(f\"  scored: {name}\", flush=True)\n",
    "for r in sorted([r for r in rows if r[\"decoder\"].startswith(\"B\")], key=lambda r: -r[\"score\"][0])[:4]:\n",
    "    rows.append(dict(scenario=r[\"scenario\"], decoder=\"C per-metro cut-points\", out=r[\"out\"], fit=fit_metro_cuts,\n",
    "                     apply=apply_metro_cuts, score=crossfit(r[\"out\"][0], fit_metro_cuts, apply_metro_cuts)))\n",
    "    print(f\"  scored per-metro: {r['scenario']}\", flush=True)\n",
    "\n",
    "table = pd.DataFrame([{\"scenario\": r[\"scenario\"], \"decoder\": r[\"decoder\"],\n",
    "                       \"CV macro F1 (test-weighted)\": round(r[\"score\"][0], 4), \"CV macro F1 (plain)\": round(r[\"score\"][1], 4)}\n",
    "                      for r in rows]).sort_values(\"CV macro F1 (test-weighted)\", ascending=False).reset_index(drop=True)\n",
    "prog.finish()\n",
    "display(table)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7008fef8",
   "metadata": {},
   "source": [
    "## 8. Submission\n",
    "`submission.csv` = the best scenario. The decoder is re-tuned on **all** out-of-fold predictions, then applied to the test predictions.\n",
    "The 2nd and 3rd best are saved too (`submission_rank2.csv`, `submission_rank3.csv`).\n",
    "**Pick by CV, not by the public leaderboard**: the public part is only a slice of the test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "176537a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "prog.start(\"write submissions\")\n",
    "ranked = sorted(rows, key=lambda r: -r[\"score\"][0])\n",
    "for i, r in enumerate(ranked[:3]):\n",
    "    o, t = r[\"out\"]\n",
    "    prm = r[\"fit\"](o, y, W_tr, MT_tr)\n",
    "    pred = r[\"apply\"](t, prm, MT_te)\n",
    "    fname = \"submission.csv\" if i == 0 else f\"submission_rank{i + 1}.csv\"\n",
    "    pd.DataFrame({\"ID\": test[\"space_id\"], \"price_tier\": pred.astype(int) + 1}).to_csv(WORK / fname, index=False)\n",
    "    print(f\"{fname:24s} \u2190 {r['scenario']} | {r['decoder']} | CV (test-weighted) {r['score'][0]:.4f}\")\n",
    "prog.finish()\n",
    "pd.read_csv(WORK / \"submission.csv\").price_tier.value_counts().sort_index()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a74d5ba9",
   "metadata": {},
   "source": [
    "## 9. Error analysis\n",
    "Where does the chosen model go wrong? The predictions below are the cross-fitted out-of-fold ones (each fold decoded with cut-points/weights tuned on the other folds)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c97ea8ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import confusion_matrix, f1_score as sk_f1\n",
    "best_r = ranked[0]\n",
    "o_, _ = best_r[\"out\"]\n",
    "cf_pred = np.empty(n_tr, int)\n",
    "for f in range(FOLDS):\n",
    "    m_ = fold_id == f\n",
    "    prm = best_r[\"fit\"](o_[~m_], y[~m_], W_tr[~m_], MT_tr[~m_])\n",
    "    cf_pred[m_] = best_r[\"apply\"](o_[m_], prm, MT_tr[m_])\n",
    "print(\"chosen:\", best_r[\"scenario\"], \"|\", best_r[\"decoder\"])\n",
    "print(\"per-tier F1:\", dict(zip([\"tier 1\", \"tier 2\", \"tier 3\", \"tier 4\"], sk_f1(y, cf_pred, average=None).round(3))))\n",
    "cm = confusion_matrix(y, cf_pred, normalize=\"true\")\n",
    "fig, ax = plt.subplots(figsize=(5, 4))\n",
    "ax.imshow(cm, cmap=\"Blues\")\n",
    "for i in range(K):\n",
    "    for j in range(K):\n",
    "        ax.text(j, i, f\"{cm[i, j]:.2f}\", ha=\"center\", va=\"center\", color=\"white\" if cm[i, j] > 0.5 else \"black\")\n",
    "ax.set_xticks(range(K), [f\"pred {k + 1}\" for k in range(K)]); ax.set_yticks(range(K), [f\"true {k + 1}\" for k in range(K)])\n",
    "ax.set_title(\"Confusion matrix (row-normalised)\"); plt.tight_layout(); plt.show()\n",
    "err = cf_pred != y\n",
    "print(f\"wrong: {err.mean():.1%} of listings | of the mistakes, {(np.abs(cf_pred - y)[err] == 1).mean():.1%} are off by just one tier\")\n",
    "per_metro = pd.DataFrame({\"metro\": train.metro.astype(str), \"y\": y, \"p\": cf_pred}).groupby(\"metro\")[[\"y\", \"p\"]].apply(\n",
    "    lambda d: pd.Series({\"listings\": len(d), \"macro F1\": round(sk_f1(d.y, d.p, average=\"macro\"), 3)}))\n",
    "display(per_metro.sort_values(\"macro F1\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e4586e1",
   "metadata": {},
   "source": [
    "**Reading this.** Look at where the mistakes sit in the confusion matrix: mostly next to the diagonal means the model ranks listings well and only the borders between tiers are blurry.\n",
    "Middle tiers usually score lowest because they border two tiers each. The per-metro table shows which markets are hardest to predict."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "baceb8c8",
   "metadata": {},
   "source": [
    "## 10. What else could be tried\n",
    "- **Pseudo-labelling**: add confident test predictions to training, so the model sees the new operators (hard to validate honestly).\n",
    "- **Seed averaging**: train each model with 3 seeds and average (small, reliable gain; costs 3\u00d7 time).\n",
    "- **Tuning**: a small random search on depth / learning rate / regularisation for XGBoost and CatBoost.\n",
    "- **Facilities as text**: TF-IDF + SVD over facility codes instead of one-hot columns.\n",
    "- **Not allowed / not used**: external data (the original public listings), AutoML, deep-learning frameworks."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "accelerator": "GPU",
  "colab": {
   "gpuType": "T4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}