{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "linclass.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# linear classifier \u2014 Python demo\n\nNumerical companion to the entry [linear classifier](https://dictionaryofml.org/terms/linclass.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nA binary classification trainset for the vineyard task: each data point is a square patch of an aerial photograph, described by two numeric features \u2014 the contrast of the patch (the standard deviation of its pixel brightness) and its greenness (the relative difference between the average green and the average red channel value, in percent) \u2014 and labelled +1 if the patch shows a vineyard and -1 otherwise. The feature vectors are drawn from a two-component Gaussian model whose means and covariance matrices are the patch statistics measured on an orthophoto of the Wachau valley in Austria (basemap.at, CC BY 4.0; 41,472 patches of 32 x 32 pixels at 0.79 m per pixel, 21 percent of them vineyard). Vine rows alternate foliage with bare soil, which makes a vineyard patch less green and more contrasted than the forest and meadow around it.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/linclass.py`](https://dictionaryofml.org/terms/linclass.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(), \"linclass.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\nlinclass.py \u2014 numerical companion to the glossary entry 'linear classifier'.\n\nPurpose\n-------\nA binary classification trainset for the vineyard task: each data point is\na square patch of an aerial photograph, described by two numeric features \u2014\nthe contrast of the patch (the standard deviation of its pixel brightness)\nand its greenness (the relative difference between the average green and\nthe average red channel value, in percent) \u2014 and labelled +1 if the patch\nshows a vineyard and -1 otherwise.  The feature vectors are drawn from a\ntwo-component Gaussian model whose means and covariance matrices are the\npatch statistics measured on an orthophoto of the Wachau valley in Austria\n(basemap.at, CC BY 4.0; 41,472 patches of 32 x 32 pixels at 0.79 m per\npixel, 21 percent of them vineyard).  Vine rows alternate foliage with\nbare soil, which makes a vineyard patch less green and more contrasted\nthan the forest and meadow around it.\n\nA linear classifier is learned from the trainset by gradient descent on\nthe average logistic loss, and the geometry of its decision boundary is\nthen checked numerically: the distance of a feature vector from the\ndecision boundary equals |h(x)| / ||w||; the hypothesis value h(x) alone\ndoes not measure that distance, since rescaling the parameters multiplies\nit by the scale factor while leaving the classifier and the distances\nunchanged; and a feature perturbation shorter than that distance never\nchanges the predicted label.  Self-contained (numpy/matplotlib only),\nfixed seed.\n\nBlocks\n------\n[B-patches] m = 240 data points, 120 per class, drawn from the two-\n            component model; the contrast is kept positive by redrawing.\n            Written to linclass_vineyard.csv and linclass_other.csv for\n            the entry's scatter plot.\n[B-fit]     The parameters (w, b) of the linear classifier h(x) = w^T x + b\n            learned by gradient descent on the average logistic loss (in\n            standardized coordinates, mapped back to the original feature\n            units).  The learned classifier beats the constant rule that\n            always answers with the majority label.  Its decision boundary\n            and its normal vector w go to linclass_boundary.csv and\n            linclass_normal.csv.\n[B-geom]    The distance of each of the 240 feature vectors from the\n            decision boundary, obtained by minimizing the Euclidean\n            distance over a fine grid of points of the boundary, agrees\n            with |h(x)| / ||w|| to four decimals.  The foot of the\n            perpendicular from one marked data point goes to\n            linclass_drop.csv.\n[B-scale]   Rescaling the parameters, (w, b) -> (c w, c b) with c = 7,\n            multiplies every hypothesis value by c but leaves every\n            prediction and every distance unchanged: h(x) by itself is\n            not a distance, h(x) / ||w|| is.\n[B-robust]  Perturbing a feature vector by a vector shorter than its\n            distance from the decision boundary never changes the\n            predicted label; a perturbation longer than that distance\n            can change it.\n[B-preview] The matplotlib preview of the scatter plot and of the share\n            of random perturbations that change a prediction, against\n            the perturbation length.\n\nOutputs\n-------\nlinclass_vineyard.csv : x1,x2 \u2014 the 120 patches labelled +1 (vineyard).\nlinclass_other.csv    : x1,x2 \u2014 the 120 patches labelled -1.\nlinclass_boundary.csv : x1,x2 \u2014 two endpoints of the decision boundary line.\nlinclass_normal.csv   : x1,x2 \u2014 base and tip of the normal-vector w arrow.\nlinclass_drop.csv     : x1,x2 \u2014 the marked data point and the foot of its\n                        perpendicular on the decision boundary.\nlinclass.png          : matplotlib preview of the figures (checking only).\n\"\"\"\n\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 = []\n\n\ndef check(name, ok):\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")\n\n\ndef sigmoid(z):\n    return 1.0 / (1.0 + np.exp(-z))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-patches]** m = 240 data points, 120 per class, drawn from the two- component model; the contrast is kept positive by redrawing. Written to linclass_vineyard.csv and linclass_other.csv for the entry's scatter plot."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# mean and covariance of (contrast, greenness) over the patches of the\n# orthophoto, by class; vineyard patches are more contrasted and less green\nMEAN = {1: np.array([18.16, 0.04]), -1: np.array([14.81, 6.39])}\nCOV = {1: np.array([[38.46, -3.04], [-3.04, 4.06]]),\n       -1: np.array([[99.44, -29.64], [-29.64, 15.84]])}\nM_PER_CLASS = 120\nrng = np.random.default_rng(4)\n\n\ndef draw(label, n):\n    \"\"\"n feature vectors of one class; the contrast is a standard deviation,\n    so a draw with a negative first feature is replaced.\"\"\"\n    out = []\n    while len(out) < n:\n        z = rng.multivariate_normal(MEAN[label], COV[label])\n        if z[0] > 0.0:\n            out.append(z)\n    return np.array(out)\n\n\nXpos, Xneg = draw(1, M_PER_CLASS), draw(-1, M_PER_CLASS)\nX = np.r_[Xpos, Xneg]\ny = np.r_[np.ones(M_PER_CLASS), -np.ones(M_PER_CLASS)]\nm = len(X)\nprint(f\"    {m} patches, {M_PER_CLASS} of each class; \"\n      f\"contrast in [{X[:, 0].min():.1f}, {X[:, 0].max():.1f}], \"\n      f\"greenness in [{X[:, 1].min():.1f}, {X[:, 1].max():.1f}]\")\ncheck(\"[B-patches] 240 data points with two features each\",\n      X.shape == (240, 2) and len(y) == 240)\ncheck(\"[B-patches] every contrast is positive\", (X[:, 0] > 0.0).all())\ncheck(\"[B-patches] the vineyard patches are less green on average\",\n      Xpos[:, 1].mean() < Xneg[:, 1].mean())\nfor fname, Xc in [(\"linclass_vineyard.csv\", Xpos), (\"linclass_other.csv\", Xneg)]:\n    with open(OUT_DIR / fname, \"w\") as f:\n        f.write(\"x1,x2\\n\")\n        for a, b_ in Xc:\n            f.write(f\"{a:.4f},{b_:.4f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-fit]** The parameters (w, b) of the linear classifier h(x) = w^T x + b learned by gradient descent on the average logistic loss (in standardized coordinates, mapped back to the original feature units). The learned classifier beats the constant rule that always answers with the majority label. Its decision boundary and its normal vector w go to linclass_boundary.csv and linclass_normal.csv."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "mu, sd = X.mean(axis=0), X.std(axis=0)\nZ = np.c_[(X - mu) / sd, np.ones(m)]      # standardized, constant feature last\n\n\ndef grad(wz):\n    return -(Z * (y * sigmoid(-y * (Z @ wz)))[:, None]).mean(axis=0)\n\n\nwz = np.zeros(3)\nfor _ in range(6000):\n    wz = wz - 0.5 * grad(wz)\nweights = wz[:2] / sd                             # w in the original units\noffset = float(wz[2] - np.sum(wz[:2] * mu / sd))  # b in the original units\nwnorm = float(np.linalg.norm(weights))\n\n\ndef hyp(pts):\n    return np.atleast_1d(pts @ weights + offset)\n\n\nacc = float(np.mean(np.sign(hyp(X)) == y))\nacc_const = float(max((y > 0).mean(), (y < 0).mean()))\nprint(f\"    w = ({weights[0]:.3f}, {weights[1]:.3f}), b = {offset:.3f}; \"\n      f\"correctly classified {acc:.3f} against {acc_const:.3f}\")\ncheck(\"[B-fit]     the learned classifier beats the constant majority rule\",\n      acc > acc_const)\ncheck(\"[B-fit]     the gradient has (almost) vanished\",\n      float(np.linalg.norm(grad(wz))) < 1e-3)\n\nX1LO, X1HI = 0.0, 48.0                            # the plotted range of x1\n\n\ndef x2_on_boundary(x1):\n    return -(weights[0] * x1 + offset) / weights[1]\n\n\nwith open(OUT_DIR / \"linclass_boundary.csv\", \"w\") as f:\n    f.write(\"x1,x2\\n\")\n    for x1 in (X1LO, X1HI):\n        f.write(f\"{x1:.4f},{x2_on_boundary(x1):.4f}\\n\")\n\n# the normal vector w, drawn from a point of the boundary into the half-space\n# that the classifier labels +1 (where h is positive)\nunit_w = weights / wnorm\nbase = np.array([42.0, x2_on_boundary(42.0)])\nARROW = 4.5                                       # arrow length in feature units\ntip = base + ARROW * unit_w\ncheck(\"[B-fit]     the arrow points into the half-space labelled +1\",\n      float(hyp(tip)[0]) > 0.0)\nwith open(OUT_DIR / \"linclass_normal.csv\", \"w\") as f:\n    f.write(\"x1,x2\\n\")\n    for pt in (base, tip):\n        f.write(f\"{pt[0]:.4f},{pt[1]:.4f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-geom]** The distance of each of the 240 feature vectors from the decision boundary, obtained by minimizing the Euclidean distance over a fine grid of points of the boundary, agrees with |h(x)| / ||w|| to four decimals. The foot of the perpendicular from one marked data point goes to linclass_drop.csv."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "dist_formula = np.abs(hyp(X)) / wnorm\n# the same distance, obtained by searching a fine grid of boundary points\nt = np.linspace(-400.0, 400.0, 400001)\nalong = np.array([-weights[1], weights[0]]) / wnorm      # unit vector, w^T . = 0\nfoot0 = -offset * weights / (weights @ weights)          # a point of the boundary\nline = foot0 + t[:, None] * along\ndist_grid = np.array([np.min(np.linalg.norm(line - x, axis=1)) for x in X])\ngap = float(np.abs(dist_grid - dist_formula).max())\nprint(f\"    largest deviation between |h(x)|/||w|| and the grid search: {gap:.2e}\")\ncheck(\"[B-geom]    |h(x)|/||w|| is the distance from the decision boundary\",\n      gap < 1e-4)\n# a vineyard patch on the sparse left flank of the cloud, as far from the\n# decision boundary as the plot can show without the drop crossing the clouds\nleft = (X[:, 0] > 5.0) & (X[:, 0] < 12.0)\nleft[M_PER_CLASS:] = False                        # vineyard patches only\nMARK = int(np.argmax(np.where(left, dist_formula, -np.inf)))\nx_mark = X[MARK]\nfoot = x_mark - float(hyp(x_mark)[0]) / (weights @ weights) * weights\nprint(f\"    marked patch ({x_mark[0]:.2f}, {x_mark[1]:.2f}): h(x) = \"\n      f\"{float(hyp(x_mark)[0]):.2f}, |h(x)|/||w|| = {dist_formula[MARK]:.2f}\")\ncheck(\"[B-geom]    the foot of the perpendicular lies on the boundary\",\n      abs(float(hyp(foot)[0])) < 1e-9)\ncheck(\"[B-geom]    the drop from the marked patch is orthogonal to the boundary\",\n      abs(float((x_mark - foot) @ along)) < 1e-9)\nwith open(OUT_DIR / \"linclass_drop.csv\", \"w\") as f:\n    f.write(\"x1,x2\\n\")\n    for pt in (x_mark, foot):\n        f.write(f\"{pt[0]:.4f},{pt[1]:.4f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-scale]** Rescaling the parameters, (w, b) -> (c w, c b) with c = 7, multiplies every hypothesis value by c but leaves every prediction and every distance unchanged: h(x) by itself is not a distance, h(x) / ||w|| is."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "C = 7.0\nh_scaled = C * hyp(X)\ndist_scaled = np.abs(h_scaled) / (C * wnorm)\ncheck(\"[B-scale]   rescaling multiplies every hypothesis value by c\",\n      np.allclose(h_scaled, C * hyp(X)) and abs(C - 1.0) > 0.5)\ncheck(\"[B-scale]   rescaling leaves every prediction unchanged\",\n      np.array_equal(np.sign(h_scaled), np.sign(hyp(X))))\ncheck(\"[B-scale]   rescaling leaves every distance unchanged\",\n      np.allclose(dist_scaled, dist_formula))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-robust]** Perturbing a feature vector by a vector shorter than its distance from the decision boundary never changes the predicted label; a perturbation longer than that distance can change it."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "n_short = n_flip_short = n_flip_long = 0\nfor idx in range(m):\n    for _ in range(20):\n        direction = rng.normal(size=2)\n        direction /= np.linalg.norm(direction)\n        short = direction * (0.99 * dist_formula[idx])\n        long_ = direction * (1.01 * dist_formula[idx])\n        n_short += 1\n        n_flip_short += int(np.sign(hyp(X[idx] + short)[0]) != np.sign(hyp(X[idx])[0]))\n        n_flip_long += int(np.sign(hyp(X[idx] + long_)[0]) != np.sign(hyp(X[idx])[0]))\nprint(f\"    {n_short} perturbations of each length: {n_flip_short} shorter ones \"\n      f\"and {n_flip_long} longer ones change the prediction\")\ncheck(\"[B-robust]  no perturbation shorter than the distance changes a prediction\",\n      n_flip_short == 0)\ncheck(\"[B-robust]  some perturbation longer than the distance does change one\",\n      n_flip_long > 0)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-preview]** The matplotlib preview of the scatter plot and of the share of random perturbations that change a prediction, against the perturbation length."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "fig, (ax, ax2) = plt.subplots(1, 2, figsize=(10.4, 4.2))\nax.plot(Xpos[:, 0], Xpos[:, 1], \"o\", color=\"0.6\", ms=5,\n        label=\"vineyard ($y = +1$)\")\nax.plot(Xneg[:, 0], Xneg[:, 1], \"^\", color=\"black\", mfc=\"none\", ms=6,\n        label=\"no vineyard ($y = -1$)\")\nax.plot([X1LO, X1HI], [x2_on_boundary(X1LO), x2_on_boundary(X1HI)], \"k-\",\n        lw=1.8, label=\"decision boundary $h(x) = 0$\")\nax.annotate(\"\", xy=tuple(tip), xytext=tuple(base),\n            arrowprops=dict(arrowstyle=\"->\", lw=1.6))\nax.text(tip[0] + 0.6, tip[1] - 0.6, \"$w$\", fontsize=12)\nax.plot([x_mark[0], foot[0]], [x_mark[1], foot[1]], \"k:\", lw=1.8)\nax.plot(*x_mark, \"o\", color=\"black\", ms=7)\nax.text(x_mark[0] + 1.0, 0.5 * (x_mark[1] + foot[1]) - 0.6,\n        \"$|h(x)|\\\\,/\\\\,\\\\|w\\\\|$\", fontsize=11)\nax.set_xlim(X1LO, X1HI)\nax.set_aspect(\"equal\")\nax.set_xlabel(\"contrast of the patch, $x_1$\")\nax.set_ylabel(\"greenness of the patch in percent, $x_2$\")\nax.set_title(\"Wachau patches: vineyard or not\")\nax.legend(frameon=False, fontsize=8, loc=\"upper right\")\nradii = np.linspace(0.4, 2.0, 17)\nflipped = []\nfor r in radii:\n    d_rand = rng.normal(size=(m, 40, 2))\n    d_rand /= np.linalg.norm(d_rand, axis=2, keepdims=True)\n    moved = X[:, None, :] + d_rand * (r * dist_formula)[:, None, None]\n    flipped.append(float(np.mean(np.sign(moved @ weights + offset)\n                                 != np.sign(hyp(X))[:, None])))\nax2.plot(radii, flipped, \"o-\", color=\"black\", ms=4)\nax2.axvline(1.0, color=\"0.6\", ls=\"--\", lw=1.4)\nax2.text(1.04, 0.28, \"perturbation length\\n$= |h(x)|\\\\,/\\\\,\\\\|w\\\\|$\", fontsize=9)\nax2.set_xlabel(\"perturbation length, in units of $|h(x)|\\\\,/\\\\,\\\\|w\\\\|$\")\nax2.set_ylabel(\"share of perturbations that change the prediction\")\nax2.set_title(\"no shorter perturbation changes a prediction\")\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"linclass.png\", dpi=110)\n\nn_ok = sum(ok for _, ok in report)\nprint(f\"\\n{n_ok}/{len(report)} checks pass\")\nprint(f\"wrote linclass_vineyard.csv, linclass_other.csv, linclass_boundary.csv, \"\n      f\"linclass_normal.csv, linclass_drop.csv, linclass.png in {OUT_DIR}\")\nif n_ok != len(report):\n    raise SystemExit(1)"
  }
 ]
}