{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "decisionregion.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# decision region \u2014 Python demo\n\nNumerical companion to the entry [decision region](https://dictionaryofml.org/terms/decisionregion.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nBacks the entry's claims numerically, on the training set and the three classifiers of the entry 'decision boundary' (pythondemos/ decisionboundary.py): the decision regions of a classifier partition the feature space, one region per label value; their shape depends on the model (two half-planes for a linear classifier, unions of axis-parallel rectangles for a decision tree, a region with a curved border for logistic regression on polynomial features); and the regions determine the classifier completely, so the prediction for a data point is the label value of the region its feature vector lies in, and the average zero-one loss is the fraction of data points lying in the region of the other label value. Self-contained (numpy and matplotlib only), no randomness.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/decisionregion.py`](https://dictionaryofml.org/terms/decisionregion.py); CC BY 4.0.\n\nBefore running a downloaded copy, scan it: `uvx bandit -q -t B102,B307,B602,B605 decisionregion.py` (B102,B307,B602,B605 are the exec, eval and shell checks of [bandit](https://bandit.readthedocs.io/)) reports nothing: the script runs no shell command and evaluates no generated code."
  },
  {
   "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(), \"decisionregion.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\ndecisionregion.py \u2014 numerical companion to the glossary entry 'decision\nregion'.\n\nPurpose\n-------\nBacks the entry's claims numerically, on the training set and the three\nclassifiers of the entry 'decision boundary' (pythondemos/\ndecisionboundary.py): the decision regions of a classifier partition the\nfeature space, one region per label value; their shape depends on the\nmodel (two half-planes for a linear classifier, unions of axis-parallel\nrectangles for a decision tree, a region with a curved border for logistic\nregression on polynomial features); and the regions determine the\nclassifier completely, so the prediction for a data point is the label\nvalue of the region its feature vector lies in, and the average zero-one\nloss is the fraction of data points lying in the region of the other label\nvalue.  Self-contained (numpy and matplotlib only), no randomness.\n\nSetup\n-----\nThe same ten data points with two features as in decisionboundary.py, whose\nlabel values y = +1 and y = -1 are separated by a parabola, and the same\nthree classifiers learned from them: logistic regression on the two\nfeatures (a linear classifier), a decision tree of depth at most four, and\nlogistic regression on all monomials of the two features up to degree two.\nThe regions are evaluated on a grid of 60 x 60 points over [-2.5, 2.5]^2.\n\nBlocks\n------\n[B-partition]  For each classifier, every grid point lies in exactly one of\n               its two decision regions, and both regions are nonempty.\n[B-shape]      The linear classifier's region of +1 is the set of grid\n               points with w^T x + w_0 >= 0, and the midpoint of two grid\n               points of one region lies in that region (half-planes,\n               convex); the decision tree's prediction is constant on each\n               rectangle cut out by its thresholds (unions of axis-parallel\n               rectangles); the lower border of the polynomial classifier's\n               region of +1 is lower at x1 = 0 than at x1 = -2 and x1 = 2\n               (not a straight line).\n[B-determined] The label value of the region containing a data point's\n               feature vector, found by the half-plane test for the linear\n               classifier and by the rectangle table for the tree, is the\n               classifier's prediction; for each classifier, the number of\n               data points lying in the region of the other label value,\n               divided by ten, is its average zero-one loss.\n\nOutputs\n-------\ndecisionregion_points.csv  : x1, x2, y of the ten data points.\ndecisionregion_regions.csv : the grid with the label predicted by the\n                             linear classifier (lin), the decision tree\n                             (tree) and the polynomial classifier (poly);\n                             60 columns per row, rows ordered by x2.\ndecisionregion.png         : matplotlib preview (checking only).\n\"\"\"\n\n# /// script\n# requires-python = \">=3.10\"\n# dependencies = [\"matplotlib==3.10.7\", \"numpy==2.3.4\"]\n# ///\n\nimport numpy as np\nimport matplotlib\n\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nfrom pathlib import Path\n\nOUT_DIR = Path(__file__).parent\n\nreport = []\n\n\ndef check(name, ok):\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")\n\n\n# the ten data points of decisionboundary.py, separated by a parabola\nX = np.array([\n    [-2.1, 1.5], [-1.2, 0.2], [0.0, 0.0], [1.2, 0.2], [2.1, 1.5],\n    [-2.1, 0.5], [-1.2, -0.5], [0.0, -1.0], [1.2, -0.5], [2.1, 0.5],\n])\ny = np.where(X[:, 1] > 0.42 * X[:, 0] ** 2 - 0.65, 1, -1)\nm = len(X)\n\n\ndef feature_map(x, degree):\n    \"\"\"Original features, or all monomials up to degree two.\"\"\"\n    if degree == 1:\n        return np.column_stack((np.ones(len(x)), x[:, 0], x[:, 1]))\n    return np.column_stack((np.ones(len(x)), x[:, 0], x[:, 1],\n                            x[:, 0] ** 2, x[:, 0] * x[:, 1], x[:, 1] ** 2))\n\n\ndef train_logistic(x, labels, degree):\n    phi = feature_map(x, degree)\n    weights = np.zeros(phi.shape[1])\n    binary_labels = (labels + 1) / 2\n    for _ in range(5000):\n        scores = np.clip(phi @ weights, -30, 30)\n        probs = 1 / (1 + np.exp(-scores))\n        weights -= 0.3 * phi.T @ (probs - binary_labels) / len(labels)\n    return weights\n\n\ndef logistic_predict(x, weights, degree):\n    return np.where(feature_map(x, degree) @ weights >= 0, 1, -1)\n\n\ndef grow_tree(x, labels, depth=0, max_depth=4, min_leaf=1):\n    \"\"\"Grow a small CART classification tree using Gini impurity.\"\"\"\n    node = {\"prediction\": 1 if np.mean(labels == 1) >= 0.5 else -1}\n    if depth == max_depth or np.all(labels == labels[0]):\n        return node\n    best = None\n    for feature in range(x.shape[1]):\n        values = np.unique(x[:, feature])\n        for cut in (values[:-1] + values[1:]) / 2:\n            left = x[:, feature] <= cut\n            if left.sum() < min_leaf or (~left).sum() < min_leaf:\n                continue\n            impurity = 0.0\n            for part in (left, ~left):\n                p = np.mean(labels[part] == 1)\n                impurity += part.mean() * 2 * p * (1 - p)\n            if best is None or impurity < best[0]:\n                best = (impurity, feature, cut, left)\n    if best is None:\n        return node\n    _, feature, cut, left = best\n    node.update(feature=feature, cut=cut,\n                left=grow_tree(x[left], labels[left], depth + 1, max_depth, min_leaf),\n                right=grow_tree(x[~left], labels[~left], depth + 1, max_depth, min_leaf))\n    return node\n\n\ndef tree_predict(x, tree):\n    predictions = np.empty(len(x), dtype=int)\n    for i, row in enumerate(x):\n        node = tree\n        while \"feature\" in node:\n            node = node[\"left\"] if row[node[\"feature\"]] <= node[\"cut\"] else node[\"right\"]\n        predictions[i] = node[\"prediction\"]\n    return predictions\n\n\ndef tree_cuts(node, acc):\n    \"\"\"The threshold values of a tree, per feature.\"\"\"\n    if \"feature\" in node:\n        acc[node[\"feature\"]].add(node[\"cut\"])\n        tree_cuts(node[\"left\"], acc); tree_cuts(node[\"right\"], acc)\n    return acc\n\n\nw_linear = train_logistic(X, y, degree=1)\nw_poly = train_logistic(X, y, degree=2)\ntree = grow_tree(X, y)\nclassifiers = [\n    (\"linear classifier\", lambda z: logistic_predict(z, w_linear, 1)),\n    (\"decision tree\", lambda z: tree_predict(z, tree)),\n    (\"polynomial classifier\", lambda z: logistic_predict(z, w_poly, 2)),\n]\n\ng = np.linspace(-2.5, 2.5, 60)\nG1, G2 = np.meshgrid(g, g)                            # row i: x2 = g[i]\nG = np.c_[G1.ravel(), G2.ravel()]\nR = {name: predict(G) for name, predict in classifiers}"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-partition]** For each classifier, every grid point lies in exactly one of its two decision regions, and both regions are nonempty."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "ok, parts = True, []\nfor name, _ in classifiers:\n    pos, neg = R[name] == 1, R[name] == -1\n    ok &= bool(np.all(pos ^ neg) and pos.any() and neg.any())\n    parts.append(f\"{name} {int(pos.sum())}/{int(neg.sum())}\")\ncheck(f\"[B-partition]  each of the {len(G)} grid points lies in exactly one \"\n      f\"of the two regions of every classifier (grid points in the region \"\n      f\"of +1 / of -1: {', '.join(parts)})\", ok)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-shape]** The linear classifier's region of +1 is the set of grid points with w^T x + w_0 >= 0, and the midpoint of two grid points of one region lies in that region (half-planes, convex); the decision tree's prediction is constant on each rectangle cut out by its thresholds (unions of axis-parallel rectangles); the lower border of the polynomial classifier's region of +1 is lower at x1 = 0 than at x1 = -2 and x1 = 2 (not a straight line)."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "r_lin = R[\"linear classifier\"]\nhalfplane = np.array_equal(feature_map(G, 1) @ w_linear >= 0, r_lin == 1)\nconvex = True\nfor label in (1, -1):\n    idx = np.flatnonzero(r_lin == label)\n    a, b = np.meshgrid(idx[::5], idx[::7])\n    mid = (G[a.ravel()] + G[b.ravel()]) / 2\n    convex &= bool(np.all(logistic_predict(mid, w_linear, 1) == label))\ncuts = tree_cuts(tree, {0: set(), 1: set()})\nedges = [np.array(sorted(cuts[j])) for j in (0, 1)]\ncell = np.searchsorted(edges[0], G[:, 0]) * (len(edges[1]) + 1) + np.searchsorted(edges[1], G[:, 1])\nrectangles = all(len(np.unique(R[\"decision tree\"][cell == c])) == 1 for c in np.unique(cell))\ncell_label = {c: int(R[\"decision tree\"][cell == c][0]) for c in np.unique(cell)}\nR_poly = R[\"polynomial classifier\"].reshape(G1.shape)\nborder = {v: g[np.argmax(R_poly[:, np.argmin(np.abs(g - v))] == 1)] for v in (-2.0, 0.0, 2.0)}\ncurved = border[0.0] < min(border[-2.0], border[2.0]) - (g[1] - g[0])\ncheck(f\"[B-shape]      linear classifier: region of +1 is the half-plane \"\n      f\"w^T x + w_0 >= 0, both regions convex; decision tree: prediction \"\n      f\"constant on each of the {len(np.unique(cell))} rectangles of its \"\n      f\"{len(edges[0])} + {len(edges[1])} thresholds; polynomial classifier: \"\n      f\"lower border of the region of +1 at x2 = {border[-2.0]:.2f}, \"\n      f\"{border[0.0]:.2f}, {border[2.0]:.2f} for x1 = -2, 0, 2\",\n      halfplane and convex and rectangles and curved)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-determined]** The label value of the region containing a data point's feature vector, found by the half-plane test for the linear classifier and by the rectangle table for the tree, is the classifier's prediction; for each classifier, the number of data points lying in the region of the other label value, divided by ten, is its average zero-one loss."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "lookup_lin = np.where(feature_map(X, 1) @ w_linear >= 0, 1, -1)\ncell_X = np.searchsorted(edges[0], X[:, 0]) * (len(edges[1]) + 1) + np.searchsorted(edges[1], X[:, 1])\nlookup_tree = np.array([cell_label[c] for c in cell_X])\nerrors = {name: int(np.sum(predict(X) != y)) for name, predict in classifiers}\ncheck(f\"[B-determined] the label value of the region containing a data \"\n      f\"point's feature vector is the classifier's prediction for all {m} \"\n      f\"data points (half-plane test, rectangle table); data points in the \"\n      f\"region of the other label value: \"\n      + \", \".join(f\"{name} {e}/{m}\" for name, e in errors.items()),\n      np.array_equal(lookup_lin, classifiers[0][1](X))\n      and np.array_equal(lookup_tree, classifiers[1][1](X))\n      and all(abs(np.mean(predict(X) != y) - errors[name] / m) < 1e-12\n              for name, predict in classifiers))\n\n# ------------------------------------------------------------------- CSV\nwith open(OUT_DIR / \"decisionregion_points.csv\", \"w\") as fh:\n    fh.write(\"x1,x2,y\\n\")\n    for a, b, c in zip(X[:, 0], X[:, 1], y):\n        fh.write(f\"{a:.3f},{b:.3f},{int(c)}\\n\")\nwith open(OUT_DIR / \"decisionregion_regions.csv\", \"w\") as fh:\n    fh.write(\"x1,x2,lin,tree,poly\\n\")\n    for a, b, c, e, f in zip(G[:, 0], G[:, 1], R[\"linear classifier\"],\n                             R[\"decision tree\"], R[\"polynomial classifier\"]):\n        fh.write(f\"{a:.3f},{b:.3f},{int(c)},{int(e)},{int(f)}\\n\")\n\n# --------------------------------------------------------------- preview\nfig, axes = plt.subplots(1, 3, figsize=(10.2, 3.5), sharex=True, sharey=True)\ntitles = {\"linear classifier\": \"linear classifier\", \"decision tree\": \"decision tree\",\n          \"polynomial classifier\": \"logistic regression\\non polynomial features\"}\nfor ax, (name, _) in zip(axes, classifiers):\n    ax.contourf(G1, G2, R[name].reshape(G1.shape), levels=[0, 1.5], colors=[\"0.85\"])\n    ax.scatter(X[y == -1, 0], X[y == -1, 1], s=22, marker=\"o\", facecolors=\"white\",\n               edgecolors=\"black\", linewidths=0.9, label=\"$y = -1$\")\n    ax.scatter(X[y == 1, 0], X[y == 1, 1], s=15, marker=\"s\", facecolors=\"black\",\n               edgecolors=\"black\", linewidths=0.5, label=\"$y = +1$\")\n    ax.set_title(f\"{titles[name]}: region of +1 shaded\", fontsize=9)\n    ax.set_xlabel(r\"feature $x_1$\"); ax.set_ylabel(r\"feature $x_2$\")\n    ax.set_aspect(\"equal\")\n    ax.set_xticks([-2, 0, 2]); ax.set_yticks([-2, 0, 2])\n    ax.spines[\"top\"].set_visible(False); ax.spines[\"right\"].set_visible(False)\naxes[0].legend(frameon=False, fontsize=7, loc=\"lower left\")\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"decisionregion.png\", dpi=110)\n\nn_ok = sum(ok for _, ok in report)\nprint(f\"\\n{n_ok}/{len(report)} checks pass\")\nprint(f\"wrote {OUT_DIR / 'decisionregion_points.csv'}, \"\n      f\"{OUT_DIR / 'decisionregion_regions.csv'}, {OUT_DIR / 'decisionregion.png'}\")\nif n_ok != len(report):\n    raise SystemExit(1)"
  }
 ]
}