{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "dataimputation.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# data imputation \u2014 Python demo\n\nNumerical companion to the entry [data imputation](https://dictionaryofml.org/terms/dataimputation.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nThe entry's construction, carried out: an aerial photograph loses a rectangular region, and the lost pixels are filled in by learning to predict the RGB values of a pixel from the RGB values of its neighbors. Every intact pixel is one data point, its feature vector holds the RGB values of the eight surrounding pixels, and its label is its own RGB triple. Self-contained (numpy/matplotlib only), fixed seed.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/dataimputation.py`](https://dictionaryofml.org/terms/dataimputation.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(), \"dataimputation.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\ndataimputation.py -- numerical companion to the glossary entry\n'data imputation'.\n\nThe entry's construction, carried out: an aerial photograph loses a\nrectangular region, and the lost pixels are filled in by learning to\npredict the RGB values of a pixel from the RGB values of its neighbors.\nEvery intact pixel is one data point, its feature vector holds the RGB\nvalues of the eight surrounding pixels, and its label is its own RGB\ntriple. Self-contained (numpy/matplotlib only), fixed seed.\n\nThe photograph is synthetic. It stands in for an aerial view of the\nRossatz area on the Danube and is built to carry what matters here: wide\nregions of nearly constant color (river, meadow, forest) separated by\nsharp boundaries, so that a neighbor is informative inside a region and\nmisleading across a boundary.\n\nBlocks\n------\n[B-scene]  The synthetic aerial view and the rectangular region whose\n           pixels are lost, with the share of pixels that went missing.\n[B-set]    Imputation read as a prediction task: every intact pixel\n           whose eight neighbors are also intact becomes a data point,\n           with a feature vector of 24 numbers and a label of 3.\n[B-fit]    Linear regression from the 24 features to the 3 label values,\n           fitted on the intact pixels, is compared on held-out intact\n           pixels against the baseline that predicts the average color.\n[B-fill]   The corrupted region is filled in by iterating the learned\n           hypothesis: the lost pixels start at the average color and are\n           re-predicted from their neighbors until the sweep stops\n           changing them, which is a fixed point of that sweep. The result\n           is scored against the pixels that were removed.\n[B-edge]   Where the error falls: the imputation is accurate inside the\n           regions of nearly constant color and worst on the pixels that\n           sit on a boundary.\n\nOutputs\n-------\ndataimputation_err.csv  : per-method average squared error, for the entry's\n                          figure.\ndataimputation.png      : preview of the scene, the gap and the filled\n                          result (checking only).\n\nData generated by pythondemos/dataimputation.py.\n\"\"\"\n\nimport numpy as np\nimport matplotlib\n\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\nfrom pathlib import Path\n\nOUT_DIR = Path(__file__).parent\nrng = np.random.default_rng(20261007)\n\nreport = []\n\n\ndef check(name, ok):\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-scene]** The synthetic aerial view and the rectangular region whose pixels are lost, with the share of pixels that went missing."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[B-scene] a synthetic aerial view and the region that is lost\")\n\nH, W = 64, 96\nrows, cols = np.mgrid[0:H, 0:W]\nscene = np.zeros((H, W, 3))\n# the Danube runs across the view, meadow above it, forest below\nriver = np.abs(rows - (26 + 6 * np.sin(cols / 14.0))) < 6\nforest = rows > 44\nmeadow = ~river & ~forest\nfor mask, colour in ((river, (0.24, 0.40, 0.55)),\n                     (meadow, (0.55, 0.68, 0.33)),\n                     (forest, (0.18, 0.34, 0.20))):\n    scene[mask] = colour\nscene += 0.025 * rng.normal(size=scene.shape)\nscene = np.clip(scene, 0.0, 1.0)\n\nGAP = (np.s_[22:38], np.s_[40:64])\nlost = np.zeros((H, W), dtype=bool)\nlost[GAP] = True\nshare = lost.mean()\nprint(f\"    view {H} x {W} pixels; the lost region covers \"\n      f\"{lost.sum()} pixels ({100 * share:.1f}%)\")\ncheck(\"[B-scene] the view carries three regions and a boundary between them\",\n      river.sum() > 0 and meadow.sum() > 0 and forest.sum() > 0)\ncheck(\"[B-scene] between 5 and 10 percent of the pixels are lost\",\n      0.05 < share < 0.10)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-set]** Imputation read as a prediction task: every intact pixel whose eight neighbors are also intact becomes a data point, with a feature vector of 24 numbers and a label of 3."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"\\n[B-set] every intact pixel with intact neighbors is a data point\")\n\nOFFSETS = [(dr, dc) for dr in (-1, 0, 1) for dc in (-1, 0, 1)\n           if (dr, dc) != (0, 0)]\n\n\ndef neighbour_features(img, r, c):\n    return np.concatenate([img[r + dr, c + dc] for dr, dc in OFFSETS])\n\n\ndef build_set(img, usable):\n    feats, labels, where = [], [], []\n    for r in range(1, H - 1):\n        for c in range(1, W - 1):\n            if not usable[r, c]:\n                continue\n            if not all(usable[r + dr, c + dc] for dr, dc in OFFSETS):\n                continue\n            feats.append(neighbour_features(img, r, c))\n            labels.append(img[r, c])\n            where.append((r, c))\n    return np.array(feats), np.array(labels), where\n\n\nintact = ~lost\nX, y, coords = build_set(scene, intact)\nprint(f\"    {len(y)} data points, {X.shape[1]} features, \"\n      f\"{y.shape[1]} label values each\")\ncheck(\"[B-set] the feature vector holds the RGB of the eight neighbors\",\n      X.shape[1] == 8 * 3)\ncheck(\"[B-set] the label is the pixel's own RGB triple\", y.shape[1] == 3)\ncheck(\"[B-set] no data point uses a lost pixel\",\n      all(intact[r, c] for r, c in coords))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-fit]** Linear regression from the 24 features to the 3 label values, fitted on the intact pixels, is compared on held-out intact pixels against the baseline that predicts the average color."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"\\n[B-fit] linear regression on the intact pixels\")\n\nperm = rng.permutation(len(y))\ncut = int(0.8 * len(y))\ntr, va = perm[:cut], perm[cut:]\n\n\ndef fit_linear(Xtr, ytr):\n    A = np.c_[Xtr, np.ones(len(Xtr))]\n    return np.linalg.lstsq(A, ytr, rcond=None)[0]\n\n\ndef predict_linear(w, Xq):\n    return np.c_[Xq, np.ones(len(Xq))] @ w\n\n\nw_hat = fit_linear(X[tr], y[tr])\npred_va = predict_linear(w_hat, X[va])\nerr_fit = float(np.mean((pred_va - y[va]) ** 2))\nmean_colour = y[tr].mean(axis=0)\nerr_base = float(np.mean((mean_colour - y[va]) ** 2))\nprint(f\"    average squared error on held-out intact pixels: \"\n      f\"{err_fit:.5f}   (predicting the average color: {err_base:.5f})\")\ncheck(\"[B-fit] the learned hypothesis beats the average-color baseline\",\n      err_fit < 0.25 * err_base)\ncheck(\"[B-fit] its error is at the scale of the noise in the view\",\n      err_fit < 0.01)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-fill]** The corrupted region is filled in by iterating the learned hypothesis: the lost pixels start at the average color and are re-predicted from their neighbors until the sweep stops changing them, which is a fixed point of that sweep. The result is scored against the pixels that were removed."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"\\n[B-fill] filling the lost region by iterating the hypothesis\")\n\n# The lost pixels start at the average color and are then re-predicted from\n# their eight current neighbors, over and over. One sweep applies the map\n# T; filling in means iterating u^(t+1) = T(u^(t)) to its fixed point.\nfilled = scene.copy()\nfilled[lost] = mean_colour\nlost_coords = [(r, c) for r in range(1, H - 1) for c in range(1, W - 1)\n               if lost[r, c]]\ncheck(\"[B-fill] no lost pixel touches the border of the view\",\n      len(lost_coords) == int(lost.sum()))\n\nTOL, MAX_SWEEPS = 1e-5, 500\nchanges = []\nfor sweep in range(MAX_SWEEPS):\n    nxt = filled.copy()\n    Xq = np.array([neighbour_features(filled, r, c) for r, c in lost_coords])\n    nxt[tuple(np.array(lost_coords).T)] = np.clip(\n        predict_linear(w_hat, Xq), 0.0, 1.0)\n    change = float(np.abs(nxt - filled).max())\n    changes.append(change)\n    filled = nxt\n    if change < TOL:\n        break\n\nerr_filled = float(np.mean((filled[lost] - scene[lost]) ** 2))\nerr_filled_base = float(np.mean((mean_colour - scene[lost]) ** 2))\nprint(f\"    {sweep + 1} sweeps to a change of {change:.2e}; \"\n      f\"average squared error on the lost pixels: {err_filled:.5f}   \"\n      f\"(average color: {err_filled_base:.5f})\")\ncheck(\"[B-fill] the sweeps reach a fixed point of that map\",\n      change < TOL and changes[-1] < changes[0])\ncheck(\"[B-fill] the filled values stay inside the range of a color\",\n      filled.min() >= 0.0 and filled.max() <= 1.0)\ncheck(\"[B-fill] filling beats predicting the average color\",\n      err_filled < 0.6 * err_filled_base)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-edge]** Where the error falls: the imputation is accurate inside the regions of nearly constant color and worst on the pixels that sit on a boundary."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"\\n[B-edge] where the error sits\")\n\nsq_err = ((filled - scene) ** 2).mean(axis=2)\nregion = np.where(river, 0, np.where(forest, 2, 1))\non_edge = np.zeros((H, W), dtype=bool)\non_edge[1:-1, 1:-1] = (\n    (region[1:-1, 1:-1][..., None]\n     != np.stack([region[1 + dr:H - 1 + dr, 1 + dc:W - 1 + dc]\n                  for dr, dc in OFFSETS], axis=-1)).any(axis=-1))\nedge_lost = lost & on_edge\nflat_lost = lost & ~on_edge\nprint(f\"    lost pixels on a boundary: {edge_lost.sum()}, error \"\n      f\"{sq_err[edge_lost].mean():.5f}\")\nprint(f\"    lost pixels inside a region: {flat_lost.sum()}, error \"\n      f\"{sq_err[flat_lost].mean():.5f}\")\ncheck(\"[B-edge] the error is larger on the boundary pixels\",\n      sq_err[edge_lost].mean() > sq_err[flat_lost].mean())\n\nwith open(OUT_DIR / \"dataimputation_err.csv\", \"w\") as fh:\n    fh.write(\"method,err\\n\")\n    fh.write(f\"average color,{err_filled_base:.5f}\\n\")\n    fh.write(f\"learned from neighbors,{err_filled:.5f}\\n\")\n    fh.write(f\"inside a region,{sq_err[flat_lost].mean():.5f}\\n\")\n    fh.write(f\"on a boundary,{sq_err[edge_lost].mean():.5f}\\n\")\n\n\n# ---------------------------------------------------------------- preview\ncorrupted = scene.copy()\ncorrupted[lost] = 1.0\n\nfig, axs = plt.subplots(1, 3, figsize=(11.0, 3.4))\nax_a, ax_b, ax_c = axs\n\nax_a.imshow(scene, interpolation=\"nearest\")\nax_a.set_xlabel(\"pixel column\")\nax_a.set_ylabel(\"pixel row\")\nax_a.set_title(\"the aerial view\", fontsize=9)\n\nax_b.imshow(corrupted, interpolation=\"nearest\")\nax_b.set_xlabel(\"pixel column\")\nax_b.set_ylabel(\"pixel row\")\nax_b.set_title(f\"{100 * share:.0f}% of the pixels lost (white)\", fontsize=9)\n\nax_c.imshow(filled, interpolation=\"nearest\")\nax_c.set_xlabel(\"pixel column\")\nax_c.set_ylabel(\"pixel row\")\nax_c.set_title(f\"filled in, squared error {err_filled:.4f}\", fontsize=9)\n\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"dataimputation.png\", dpi=110)\n\nn_ok = sum(ok for _, ok in report)\nprint(f\"\\n{n_ok}/{len(report)} checks pass\")\nprint(\"wrote dataimputation_err.csv, dataimputation.png\")\nif n_ok != len(report):\n    raise SystemExit(1)"
  }
 ]
}