{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "featureimportance.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# feature importance \u2014 Python demo\n\nNumerical companion to the entry [feature importance](https://dictionaryofml.org/terms/featureimportance.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nNumerical companion to the glossary entry 'featureimportance'. The GeoSphere Austria weather station Krems (station id 3805, 48.42 N, 15.62 E) records six daily observations: maximum and minimum air temperature, mean air pressure, mean relative humidity, sunshine duration and precipitation. This script downloads the records for all 366 days of 2024 (dataset klima-v2-1d) and predicts the maximum daytime temperature of the next day from the twelve observations of the current and the previous day. A linear hypothesis learned by linear regression on the standardized features reaches an average squared error loss of 8.78 (squared degrees C) against a label variance of 84.59. Two importance scores are then computed for each of the twelve features: the magnitude of its learned weight, and its permutation importance (the increase of the average squared error loss when the feature's values are randomly permuted across the data points). Both scores agree: today's maximum temperature carries almost all of the credit, and the permutation importance matches its theoretical value of twice the squared weight. Today's minimum temperature obtains a weight magnitude of 0.02 and a permutation importance of 0.00: the predictions do not rest on it.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/featureimportance.py`](https://dictionaryofml.org/terms/featureimportance.py); CC BY 4.0."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# Notebook shim: the script resolves output paths relative to __file__,\n# which a notebook kernel does not define; everything lands in the\n# working directory instead.\nimport os\n__file__ = os.path.join(os.getcwd(), \"featureimportance.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"Feature importance for a next-day temperature forecast at Krems an der\nDonau: which of today's and the previous day's weather observations does\na learned linear hypothesis actually use?\n\nPurpose\n-------\nNumerical companion to the glossary entry 'featureimportance'.  The\nGeoSphere Austria weather station Krems (station id 3805, 48.42 N,\n15.62 E) records six daily observations: maximum and minimum air\ntemperature, mean air pressure, mean relative humidity, sunshine\nduration and precipitation.  This script downloads the records for all\n366 days of 2024 (dataset klima-v2-1d) and predicts the maximum\ndaytime temperature of the next day from the twelve observations of\nthe current and the previous day.  A linear hypothesis learned by\nlinear regression on the standardized features reaches an average squared\nerror loss of 8.78 (squared degrees C) against a label variance of\n84.59.  Two importance scores are then computed for each of the twelve\nfeatures: the magnitude of its learned weight, and its permutation\nimportance (the increase of the average squared error loss when the\nfeature's values are randomly permuted across the data points).  Both\nscores agree: today's maximum temperature carries almost all of the\ncredit, and the permutation importance matches its theoretical value\nof twice the squared weight.  Today's minimum temperature obtains a\nweight magnitude of 0.02 and a permutation importance of 0.00: the\npredictions do not rest on it.\n\nDeterministic: the linear hypothesis is the closed-form solution of\nlinear regression; the permutations use a fixed seed.  Self-contained: numpy +\nmatplotlib only (stdlib urllib for the download).\n\nBlocks\n------\n[B-fetch]   Download the six daily observations at Krems for 2024 from\n            the GeoSphere data hub; write the 366 records to\n            featureimportance_weather.csv and pin them to the archive.\n[B-model]   Build 364 data points (features: the twelve observations of\n            the current and the previous day, standardized; label: the\n            maximum daytime temperature of the next day); learn a\n            linear hypothesis by linear regression; write the actual and\n            predicted values for June 2024 to featureimportance_stem.csv.\n[B-weights] Score each feature by the magnitude of its learned weight;\n            today's maximum temperature dominates with 8.60, today's\n            minimum temperature is last with 0.02.\n[B-permute] Score each feature by its permutation importance (20\n            permutation rounds, fixed seed); check that it matches\n            twice the squared weight for the three largest weights;\n            write both scores to featureimportance_importance.csv.\n\nOutputs\n-------\nfeatureimportance_weather.csv    : date + six observations, 366 days\nfeatureimportance_stem.csv       : day, actual, predicted (June 2024)\nfeatureimportance_importance.csv : idx, label, wabs, pisqrt (ascending\n                                   by weight magnitude, so the largest\n                                   score sits on top of an xbar chart)\nfeatureimportance.png            : preview (checking only) -- the June\n                                   stem plot and both importance scores\n\"\"\"\n\nimport json\nimport urllib.request\nfrom pathlib import Path\n\nimport numpy as np\nimport matplotlib\n\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nOUT_DIR = Path(__file__).parent\n\nreport = []                         # collects (check name, pass/fail) pairs\n\n\ndef check(name, ok):                # records and prints one verification\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-fetch]** Download the six daily observations at Krems for 2024 from the GeoSphere data hub; write the 366 records to featureimportance_weather.csv and pin them to the archive."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "PARAMS = [\"tlmax\", \"tlmin\", \"p_mittel\", \"rf_mittel\", \"so_h\", \"rr\"]\nURL = (\"https://dataset.api.hub.geosphere.at/v1/station/historical/\"\n       f\"klima-v2-1d?parameters={','.join(PARAMS)}&station_ids=3805\"\n       \"&start=2024-01-01&end=2024-12-31\")\nwith urllib.request.urlopen(URL, timeout=120) as resp:\n    payload = json.load(resp)\ncols = payload[\"features\"][0][\"properties\"][\"parameters\"]\ndays = [stamp[:10] for stamp in payload[\"timestamps\"]]\nraw = np.array([cols[p][\"data\"] for p in PARAMS], dtype=float).T\nwith open(OUT_DIR / \"featureimportance_weather.csv\", \"w\") as f:\n    f.write(\"date,\" + \",\".join(PARAMS) + \"\\n\")\n    for day, row in zip(days, raw):\n        f.write(day + \",\" + \",\".join(f\"{v:g}\" for v in row) + \"\\n\")\ncheck(\"[B-fetch] 366 days downloaded for 2024\", len(days) == 366)\ncheck(\"[B-fetch] no observation is missing\", not np.isnan(raw).any())\ncheck(\"[B-fetch] the record matches the archive (Jan 1: 9.6, 2.6, \"\n      \"988.3, 67, 4.0, 0.0)\",\n      raw[0].tolist() == [9.6, 2.6, 988.3, 67.0, 4.0, 0.0])\nraw[:, 5] = np.maximum(raw[:, 5], 0.0)   # the archive codes trace\n# precipitation as -1.0 mm; a negative rainfall amount is a code, not a\n# measurement, so it is set to zero"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-model]** Build 364 data points (features: the twelve observations of the current and the previous day, standardized; label: the maximum daytime temperature of the next day); learn a linear hypothesis by linear regression; write the actual and predicted values for June 2024 to featureimportance_stem.csv."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "OBS = [\"max temp\", \"min temp\", \"pressure\", \"humidity\", \"sunshine\", \"precip\"]\nlabels = ([f\"{o} (today)\" for o in OBS]\n          + [f\"{o} (prev day)\" for o in OBS])\nfeats = np.array([np.concatenate([raw[t], raw[t - 1]])\n                  for t in range(1, len(raw) - 1)])\ny = raw[2:, 0]                       # tlmax of the following day\nfeats = (feats - feats.mean(0)) / feats.std(0)\nA = np.hstack([feats, np.ones((len(feats), 1))])\nw, *_ = np.linalg.lstsq(A, y, rcond=None)\npred = A @ w\nmse = float(np.mean((y - pred) ** 2))\nvar = float(np.mean((y - y.mean()) ** 2))\nprint(f\"  average squared error loss {mse:.2f}, label variance {var:.2f}\")\ncheck(\"[B-model] 364 data points with 12 features each\",\n      feats.shape == (364, 12))\ncheck(\"[B-model] the learned hypothesis reaches average loss 8.78\",\n      round(mse, 2) == 8.78)\ncheck(\"[B-model] the label variance is 84.59\", round(var, 2) == 84.59)\n\njune = [i for i, t in enumerate(range(1, len(raw) - 1))\n        if days[t + 1].startswith(\"2024-06\")]\nwith open(OUT_DIR / \"featureimportance_stem.csv\", \"w\") as f:\n    f.write(\"day,actual,predicted\\n\")\n    for n, i in enumerate(june, start=1):\n        f.write(f\"{n},{y[i]:.1f},{pred[i]:.1f}\\n\")\ncheck(\"[B-model] the stem plot covers the 30 days of June 2024\",\n      len(june) == 30)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-weights]** Score each feature by the magnitude of its learned weight; today's maximum temperature dominates with 8.60, today's minimum temperature is last with 0.02."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "wabs = np.abs(w[:12])\norder = np.argsort(wabs)             # ascending: largest ends up on top\nfor i in order[::-1]:\n    print(f\"  |weight| {labels[i]:22s} {wabs[i]:5.2f}\")\ncheck(\"[B-weights] the largest weight magnitude is today's max temp, 8.60\",\n      labels[int(order[-1])] == \"max temp (today)\"\n      and round(float(wabs[order[-1]]), 2) == 8.60)\ncheck(\"[B-weights] the smallest weight magnitude is today's min temp, 0.02\",\n      labels[int(order[0])] == \"min temp (today)\"\n      and round(float(wabs[order[0]]), 2) == 0.02)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-permute]** Score each feature by its permutation importance (20 permutation rounds, fixed seed); check that it matches twice the squared weight for the three largest weights; write both scores to featureimportance_importance.csv."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "rng = np.random.default_rng(0)\nperm = np.zeros(12)\nfor j in range(12):\n    inc = []\n    for _ in range(20):\n        Ap = A.copy()\n        Ap[:, j] = rng.permutation(Ap[:, j])\n        inc.append(float(np.mean((y - Ap @ w) ** 2)) - mse)\n    perm[j] = np.mean(inc)\npisqrt = np.sqrt(np.maximum(perm, 0.0) / 2.0)\nfor i in order[::-1]:\n    print(f\"  perm. importance {labels[i]:22s} {perm[i]:7.2f}\")\ncheck(\"[B-permute] the largest permutation importance is today's max \"\n      \"temp, 147.21\", round(float(perm[int(order[-1])]), 2) == 147.21)\ncheck(\"[B-permute] today's min temp has permutation importance 0.00\",\n      round(float(perm[int(order[0])]), 2) == 0.00)\ncheck(\"[B-permute] the importance matches twice the squared weight \"\n      \"(within 0.1 after sqrt) for the three largest weights\",\n      all(abs(pisqrt[i] - wabs[i]) < 0.1 for i in order[-3:]))\nwith open(OUT_DIR / \"featureimportance_importance.csv\", \"w\") as f:\n    f.write(\"idx,label,wabs,pisqrt\\n\")\n    for n, i in enumerate(order):\n        f.write(f\"{n},{labels[i]},{wabs[i]:.2f},{pisqrt[i]:.2f}\\n\")\n\n# ---- preview figure (checking only; the entry reads the CSVs)\nfig, (ax1, ax2) = plt.subplots(2, 1, figsize=(7.0, 7.6))\nd = np.arange(1, 31)\nya = np.array([y[i] for i in june])\nyp = np.array([pred[i] for i in june])\nax1.vlines(d, 0, ya, color=\"0.6\", linewidth=1)\nax1.plot(d, ya, \"o\", color=\"0.3\", label=\"actual\")\nax1.plot(d, yp, \"x\", color=\"black\", markersize=7, label=\"predicted\")\nax1.set_xlabel(\"day in June 2024\")\nax1.set_ylabel(\"max daytime temp (deg C)\")\nax1.set_title(\"Krems 2024: actual vs. predicted next-day max temperature\")\nax1.legend(frameon=False)\npos = np.arange(12)\nax2.barh(pos + 0.19, wabs[order], height=0.38, color=\"0.6\",\n         label=\"weight magnitude\")\nax2.barh(pos - 0.19, pisqrt[order], height=0.38, color=\"white\",\n         edgecolor=\"black\", hatch=\"///\",\n         label=\"(perm. importance / 2)$^{1/2}$\")\nax2.set_yticks(pos)\nax2.set_yticklabels([labels[i] for i in order], fontsize=8)\nax2.set_xlabel(\"importance score (deg C per standard deviation)\")\nax2.set_ylabel(\"feature\")\nax2.set_title(\"Two importance scores for the twelve features\")\nax2.legend(frameon=False)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"featureimportance.png\", dpi=110)\n\nfailed = [name for name, ok in report if not ok]\nprint(f\"{len(report) - len(failed)}/{len(report)} checks passed\"\n      + (f\"; FAILED: {failed}\" if failed else \"\"))"
  }
 ]
}