{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "66728f37",
   "metadata": {},
   "source": [
    "# Ch 00 — Math & Python Prereqs (notebook)\n",
    "\n",
    "`[start of curriculum]` · **this notebook** · `[01 the-ml-landscape →]`\n",
    "\n",
    "Runs top-to-bottom in ~1 min on free Colab CPU. Last verified 2026-06-11.\n",
    "\n",
    "**What you'll build**\n",
    "- A finite-difference derivative checker that disagrees with the analytic gradient until you fix the step size, the way Karpathy opens micrograd.\n",
    "- A fully-connected layer's forward pass `Y = X @ W.T + b` from raw shapes, plus the broadcast rules that make it one line.\n",
    "- The least-squares gradient `X.T @ (X w - y)` derived on paper, then confirmed against `torch.autograd` to five decimals.\n",
    "- `softmax`, `cross_entropy`, and `entropy` from scratch, checked against `torch.nn.functional`, including the overflow bug that makes a naive softmax return `nan`.\n",
    "- A from-scratch gradient descent that recovers `y = 3x + 5` from noisy data, the smallest end-to-end ML algorithm.\n",
    "\n",
    "**How every notebook in this curriculum works (the house protocol).** This is the first notebook, so it states the rules once. Code cells that contain a `# TODO` are yours to fill in. The last lines of that same cell are self-checks: run the cell to grade yourself. `[ ok ]` passed, `[FAIL]` prints what went wrong and how to fix it, `[ -- ]` means you have not attempted it yet. Below each exercise is a hint ladder in fold-out `<details>` blocks (open only as many rungs as you need) and then a folded solution cell that redefines the function correctly and re-checks it. Because the solution cell runs, the notebook completes top-to-bottom even if you fill in nothing. The checks always grade YOUR code in the live namespace, never a hidden answer key. There are no hidden tests and nothing is graded by an LLM. Section 1 below walks the protocol in detail.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8321b36e",
   "metadata": {},
   "source": [
    "## Before you start\n",
    "\n",
    "Three predictions to set intention. Answer before you run anything; the rest of the notebook is the reveal.\n",
    "\n",
    "1. You have a NumPy array of shape `(3, 4)` and you add an array of shape `(4,)`. Does it work, and what shape comes out? <details><summary>Answer</summary>It works. The `(4,)` array is stretched across the 3 rows to `(3, 4)`. Broadcasting aligns shapes from the right; a trailing axis of size 4 matches a length-4 vector. Part 1 makes this mechanical.</details>\n",
    "2. `f(x) = x**4 - 3*x**2 + 2`. What is `f'(x)`, and roughly where are its minima? <details><summary>Answer</summary>`f'(x) = 4*x**3 - 6*x`, zero at `x = 0` and `x = ±sqrt(1.5) ≈ ±1.22`. The two outer roots are minima, the middle one a local max. Part 3 checks this with a numerical derivative.</details>\n",
    "3. A fair coin: what is the entropy in nats, and what would it be if the coin always landed heads? <details><summary>Answer</summary>Fair coin: `-0.5*ln(0.5) - 0.5*ln(0.5) = ln 2 ≈ 0.693` nats. A certain outcome has entropy 0: no surprise. Part 5 computes both.</details>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bc5dd70",
   "metadata": {},
   "source": [
    "## Setup\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c2188b0d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:01.448112Z",
     "iopub.status.busy": "2026-06-10T18:44:01.448012Z",
     "iopub.status.idle": "2026-06-10T18:44:02.306585Z",
     "shell.execute_reply": "2026-06-10T18:44:02.306247Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "numpy 2.2.6 · torch 2.12.0+cpu · device cpu\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import torch\n",
    "print(f\"numpy {np.__version__} · torch {torch.__version__} · device cpu\")\n",
    "if np.__version__ < \"2.0\":\n",
    "    print(\"WARN: this notebook is written for NumPy 2.x; older versions may differ slightly\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "43daebff",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.307771Z",
     "iopub.status.busy": "2026-06-10T18:44:02.307648Z",
     "iopub.status.idle": "2026-06-10T18:44:02.315345Z",
     "shell.execute_reply": "2026-06-10T18:44:02.314928Z"
    }
   },
   "outputs": [],
   "source": [
    "import os, random\n",
    "SEED = 0\n",
    "FAST = bool(os.environ.get('NB_FAST'))   # CI smoke mode: ~10x fewer steps, same code paths\n",
    "STEPS = 500 if FAST else 5000            # GD steps; expected loss documented in the experiment log\n",
    "rng = np.random.default_rng(SEED)        # the one RNG we pass around\n",
    "torch.manual_seed(SEED); random.seed(SEED)\n",
    "\n",
    "# ── house self-check harness (identical across all chapter notebooks) ──\n",
    "import numpy as _np\n",
    "\n",
    "def check(label, test_fn, required=False):\n",
    "    \"\"\"Run one self-check. test_fn raises AssertionError (with a teaching\n",
    "    message) on failure, NotImplementedError if the stub is unfilled.\n",
    "    required=True is used only in solution cells; it is what CI grades.\"\"\"\n",
    "    try:\n",
    "        test_fn()\n",
    "    except NotImplementedError:\n",
    "        if required:\n",
    "            raise AssertionError(f\"{label}: reference solution incomplete\")\n",
    "        print(f\"[ -- ] {label}: not attempted yet — fill in the TODO above, then re-run.\")\n",
    "        return False\n",
    "    except AssertionError as e:\n",
    "        if required:\n",
    "            raise\n",
    "        print(f\"[FAIL] {label}: {e}\")\n",
    "        return False\n",
    "    print(f\"[ ok ] {label}\")\n",
    "    return True\n",
    "\n",
    "def attempted(*vals):\n",
    "    \"\"\"Treat None placeholders as 'not attempted'.\"\"\"\n",
    "    if any(v is None for v in vals):\n",
    "        raise NotImplementedError\n",
    "\n",
    "def check_shape(x, want):\n",
    "    assert tuple(x.shape) == tuple(want), \\\n",
    "        f\"shape {tuple(x.shape)}, expected {tuple(want)} — check your reshape/transpose order\"\n",
    "\n",
    "def check_close(got, want, atol=1e-5, rtol=1e-4, msg=\"\"):\n",
    "    g, w = _np.asarray(got, dtype=float), _np.asarray(want, dtype=float)\n",
    "    assert g.shape == w.shape, f\"shape {g.shape} vs expected {w.shape}. {msg}\"\n",
    "    bad = ~_np.isclose(g, w, atol=atol, rtol=rtol)\n",
    "    assert not bad.any(), \\\n",
    "        f\"{bad.mean():.2%} of values wrong (max diff {abs(g - w).max():.3g}). {msg}\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfa89dfd",
   "metadata": {},
   "source": [
    "> **Note:** seeds make this notebook's printed numbers reproduce on CPU. Library versions and BLAS threading can shift the last digit or two; quoted numbers hold for the pinned environment. If a value reads 0.6931 here and 0.6932 on your machine, you did nothing wrong. We re-seed at the top of any cell that draws random numbers, so a mid-notebook re-run reproduces the same output.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd940ea6",
   "metadata": {},
   "source": [
    "## The map\n",
    "\n",
    "> **Part 1 — Arrays, shapes, broadcasting.** The `np.ndarray` and the two ideas that trip everyone up: the `axis` argument and broadcasting. You build a fully-connected layer's forward pass from shapes alone.\n",
    "> **Part 2 — Linear algebra, the four operations.** Dot product, matrix-vector, matrix-matrix, norms. You implement matmul from scratch and prove it matches `@`.\n",
    "> **Part 3 — Calculus and the gradient.** A numerical derivative the Karpathy way, the chain rule, then the least-squares gradient verified against `torch.autograd`.\n",
    "> **Part 4 — Probability and maximum likelihood.** Expectation, variance, Bayes, and why negative log-likelihood is the loss behind almost everything.\n",
    "> **Part 5 — Information theory.** Entropy, cross-entropy, KL. You build `softmax` and `cross_entropy` from scratch, break a naive softmax on purpose, and fix it.\n",
    "> **Part 6 — Putting it together.** Gradient descent from scratch recovers `y = 3x + 5`. The smallest complete ML algorithm, with an experiment log.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6730615e",
   "metadata": {},
   "source": [
    "## Part 1 — Arrays, shapes, broadcasting\n",
    "\n",
    "> **Objectives.** Construct and inspect `np.ndarray`s; read a shape and predict the result of an operation; understand `axis` and broadcasting well enough to build a fully-connected layer in one line.\n",
    "\n",
    "The one object every Python ML library is built on is `np.ndarray`: a multi-dimensional, homogeneous, typed buffer with a shape. Everything else is operations on shapes. We seed the RNG first so the printed numbers reproduce.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bb0e4fdb",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.316405Z",
     "iopub.status.busy": "2026-06-10T18:44:02.316287Z",
     "iopub.status.idle": "2026-06-10T18:44:02.319205Z",
     "shell.execute_reply": "2026-06-10T18:44:02.318692Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a.shape=(4,) dtype=int64\n",
      "r.shape=(3, 4) ndim=2 size=12\n"
     ]
    }
   ],
   "source": [
    "rng = np.random.default_rng(SEED)   # re-seed so this cell reproduces in isolation\n",
    "a = np.array([1, 2, 3, 4])          # shape (4,), dtype int64\n",
    "b = np.zeros((3, 4))                # shape (3, 4), dtype float64\n",
    "r = rng.standard_normal((3, 4))     # standard normal, shape (3, 4)\n",
    "print(f\"a.shape={a.shape} dtype={a.dtype}\")\n",
    "print(f\"r.shape={r.shape} ndim={r.ndim} size={r.size}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc73fbd1",
   "metadata": {},
   "source": [
    "> **Notice that** a 1D array has shape `(4,)` with a trailing comma, while `r` has shape `(3, 4)`. The number of integers in the shape tuple is `ndim`; their product is `size`. Almost every bug in array code is a shape you did not expect, so the fastest debugging habit in ML is to print shapes everywhere.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bb844a6",
   "metadata": {},
   "source": [
    "### Reductions and the `axis` argument\n",
    "\n",
    "A reduction collapses one axis. `r.sum(axis=0)` collapses axis 0 (the rows), leaving one value per column, so the result has shape `(4,)`. `r.sum(axis=1)` collapses the columns, leaving one per row, shape `(3,)`. Get this backwards and your loss curve looks fine while the gradient is for the wrong quantity.\n",
    "\n",
    "> **Predict:** for `r` of shape `(3, 4)`, what shape is `r.sum(axis=0)`? What about `r.sum(axis=1)`? Run the cell to check.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "546bd258",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.320068Z",
     "iopub.status.busy": "2026-06-10T18:44:02.319996Z",
     "iopub.status.idle": "2026-06-10T18:44:02.322040Z",
     "shell.execute_reply": "2026-06-10T18:44:02.321791Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "full sum (scalar): () -> 0.2648\n",
      "axis=0 (per column): (4,)\n",
      "axis=1 (per row):    (3,)\n",
      "axis=1, keepdims:    (3, 1)\n"
     ]
    }
   ],
   "source": [
    "print(\"full sum (scalar):\", r.sum().shape, \"->\", round(float(r.sum()), 4))\n",
    "print(\"axis=0 (per column):\", r.sum(axis=0).shape)   # collapses 3 rows -> (4,)\n",
    "print(\"axis=1 (per row):   \", r.sum(axis=1).shape)   # collapses 4 cols -> (3,)\n",
    "print(\"axis=1, keepdims:   \", r.sum(axis=1, keepdims=True).shape)  # (3, 1), kept for broadcasting"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da98a122",
   "metadata": {},
   "source": [
    "> **Common confusion:** `r.sum(axis=1)` has shape `(3,)`, but `r.sum(axis=1, keepdims=True)` has shape `(3, 1)`. The second is what you want when you will broadcast the result back against `r`, for example to subtract each row's mean. A `(3,)` vector broadcasts against the columns, not the rows, which is almost never what you meant.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cc65cc7",
   "metadata": {},
   "source": [
    "### Broadcasting\n",
    "\n",
    "When you combine two arrays of different shapes, NumPy aligns their shapes from the right and stretches any axis of size 1 (or absent) to match. `(3, 4) + (4,)` stretches `(4,)` across the 3 rows. `(3, 1) + (1, 4)` stretches both, giving `(3, 4)`: an outer sum. This is the operation that makes vectorized code short, and the operation that produces silent shape bugs when you forget that `(3,)` and `(3, 1)` are different.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "15f17a69",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.322690Z",
     "iopub.status.busy": "2026-06-10T18:44:02.322625Z",
     "iopub.status.idle": "2026-06-10T18:44:02.325284Z",
     "shell.execute_reply": "2026-06-10T18:44:02.324952Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shape (3, 2)\n",
      "[[11. 21.]\n",
      " [12. 22.]\n",
      " [13. 23.]]\n",
      "[ ok ] broadcasting (3,1)+(2,) -> (3,2)\n"
     ]
    }
   ],
   "source": [
    "col = np.array([[1.0], [2.0], [3.0]])   # shape (3, 1)\n",
    "row = np.array([10.0, 20.0])             # shape (2,)\n",
    "outer = col + row                        # (3, 1) + (2,) -> (3, 2)\n",
    "print(\"shape\", outer.shape)\n",
    "print(outer)\n",
    "assert outer[0, 0] == 11.0 and outer[2, 1] == 23.0, \\\n",
    "    \"broadcast result wrong; (3,1)+(2,) should add each row value to each col value\"\n",
    "print(\"[ ok ] broadcasting (3,1)+(2,) -> (3,2)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c59b83c4",
   "metadata": {},
   "source": [
    "> **Interpretation.** `col` has a real column axis of length 3 and a length-1 axis that gets stretched to 2; `row` has no leading axis, so it gains one of length 3. The result is every pairwise sum. That `assert` directly under the claim is the house style: a claim that can be a check should be a check.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "351f7ee5",
   "metadata": {},
   "source": [
    "### Exercise 0.1 — Center each row with broadcasting\n",
    "`Difficulty 1/5 · ~5 min`\n",
    "\n",
    "Fill in `center_rows(X)` so it subtracts each row's mean from that row, with no Python loop. The catch is the `keepdims` from two cells ago. The check verifies every row now sums to (almost) zero and that you did not change the shape.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e1d2017b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.325990Z",
     "iopub.status.busy": "2026-06-10T18:44:02.325928Z",
     "iopub.status.idle": "2026-06-10T18:44:02.329775Z",
     "shell.execute_reply": "2026-06-10T18:44:02.329442Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ -- ] 0.1 center_rows: not attempted yet — fill in the TODO above, then re-run.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def center_rows(X):\n",
    "    \"\"\"Subtract each row's mean from that row. X: (n, d) -> (n, d).\"\"\"\n",
    "    X = np.asarray(X, dtype=float)\n",
    "    # TODO 1: per-row mean, shape (n, 1) so it broadcasts back against X (which keyword?)\n",
    "    row_means = None\n",
    "    attempted(row_means)\n",
    "    # TODO 2: subtract and return\n",
    "    return X - row_means\n",
    "\n",
    "def _centered_rows_sum_to_zero():\n",
    "    Z = center_rows(rng.standard_normal((5, 7)))\n",
    "    check_shape(Z, (5, 7))\n",
    "    check_close(Z.sum(axis=1), np.zeros(5), atol=1e-9,\n",
    "                msg=\"each row should sum to 0 after centering; did you use keepdims=True?\")\n",
    "\n",
    "check(\"0.1 center_rows\", _centered_rows_sum_to_zero)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a7fdc83",
   "metadata": {},
   "source": [
    "<details><summary>Hint 1 (conceptual)</summary>You need the mean along `axis=1`. To subtract it row-wise, the mean must keep its column axis so it broadcasts back over `X`. That is the `keepdims` argument.</details>\n",
    "\n",
    "<details><summary>Hint 2 (pseudocode)</summary>\n",
    "\n",
    "```python\n",
    "row_means = X.mean(axis=..., keepdims=...)   # shape (n, 1)\n",
    "return X - row_means\n",
    "```\n",
    "</details>\n",
    "\n",
    "<details><summary>Help — \"operands could not be broadcast together\" or rows don't sum to zero</summary>You probably wrote `X.mean(axis=1)` without `keepdims=True`. That gives shape `(n,)`, which broadcasts against the columns, subtracting the wrong thing (or erroring). Print `X.mean(axis=1).shape` and `X.mean(axis=1, keepdims=True).shape` to see the difference.</details>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c1654048",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.330431Z",
     "iopub.status.busy": "2026-06-10T18:44:02.330364Z",
     "iopub.status.idle": "2026-06-10T18:44:02.333053Z",
     "shell.execute_reply": "2026-06-10T18:44:02.332807Z"
    },
    "collapsed": true,
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ ok ] 0.1 center_rows\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#@title Solution { display-mode: \"form\" }\n",
    "# solution: redefines center_rows; the check below re-verifies the reference.\n",
    "def center_rows(X):\n",
    "    X = np.asarray(X, dtype=float)\n",
    "    row_means = X.mean(axis=1, keepdims=True)   # (n, 1) so it broadcasts back over X\n",
    "    return X - row_means\n",
    "\n",
    "check(\"0.1 center_rows\", _centered_rows_sum_to_zero, required=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40153525",
   "metadata": {},
   "source": [
    "### The fully-connected layer is one broadcast\n",
    "\n",
    "A linear layer maps a batch of inputs to a batch of outputs:\n",
    "$$\\mathbf{Y} = X W^\\top + \\mathbf{b}$$\n",
    "where $X \\in \\mathbb{R}^{B \\times d_{\\text{in}}}$, $W \\in \\mathbb{R}^{d_{\\text{out}} \\times d_{\\text{in}}}$, $\\mathbf{b} \\in \\mathbb{R}^{d_{\\text{out}}}$, and $\\mathbf{Y} \\in \\mathbb{R}^{B \\times d_{\\text{out}}}$. The matmul handles the per-example dot products; the bias is a length-$d_{\\text{out}}$ vector broadcast across all $B$ rows. Same `XW^T` you will write a hundred times in the deep-learning chapters.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "34c71d03",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.333688Z",
     "iopub.status.busy": "2026-06-10T18:44:02.333623Z",
     "iopub.status.idle": "2026-06-10T18:44:02.335764Z",
     "shell.execute_reply": "2026-06-10T18:44:02.335466Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Y.shape (32, 50)\n",
      "[ ok ] fully-connected forward (32,100) -> (32,50)\n"
     ]
    }
   ],
   "source": [
    "B, d_in, d_out = 32, 100, 50\n",
    "X = rng.standard_normal((B, d_in))\n",
    "W = rng.standard_normal((d_out, d_in))   # rows are output units\n",
    "bias = rng.standard_normal(d_out)        # one bias per output unit\n",
    "Y = X @ W.T + bias                       # (32,100) @ (100,50) + (50,) -> (32, 50)\n",
    "print(\"Y.shape\", Y.shape)\n",
    "assert Y.shape == (B, d_out), \"fully-connected output should be (batch, d_out)\"\n",
    "print(\"[ ok ] fully-connected forward (32,100) -> (32,50)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "623c634d",
   "metadata": {},
   "source": [
    "> **Key takeaways.** A shape is a tuple; `ndim` counts its entries. A reduction collapses an `axis`; `keepdims=True` preserves it so the result still broadcasts. Broadcasting aligns shapes from the right and stretches size-1 axes. The whole forward pass of a dense layer is one matmul plus a broadcast add. Print shapes when confused.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d5033b7",
   "metadata": {},
   "source": [
    "## Part 2 — Linear algebra, the four operations\n",
    "\n",
    "> **Objectives.** Use the four operations ML actually needs (dot product, matrix-vector, matrix-matrix, norms), and prove your hand-written matmul matches NumPy's `@` with an agreement check.\n",
    "\n",
    "The 800-page textbook covers eight semesters. Machine learning uses about four operations. First, the inner (dot) product of two vectors, $\\mathbf{x}^\\top \\mathbf{y} = \\sum_i x_i y_i$. It is a neuron's pre-activation $\\mathbf{w}^\\top \\mathbf{x}$, an attention score $\\mathbf{q}^\\top \\mathbf{k}$, the numerator of cosine similarity.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9174598d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.336464Z",
     "iopub.status.busy": "2026-06-10T18:44:02.336403Z",
     "iopub.status.idle": "2026-06-10T18:44:02.338528Z",
     "shell.execute_reply": "2026-06-10T18:44:02.338269Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x @ y        = 32.0\n",
      "np.dot(x, y) = 32.0\n",
      "(x*y).sum()  = 32.0\n",
      "[ ok ] dot product = 32\n"
     ]
    }
   ],
   "source": [
    "x = np.array([1., 2., 3.])\n",
    "y = np.array([4., 5., 6.])\n",
    "print(\"x @ y        =\", x @ y)          # 1*4 + 2*5 + 3*6 = 32\n",
    "print(\"np.dot(x, y) =\", np.dot(x, y))   # same\n",
    "print(\"(x*y).sum()  =\", (x * y).sum())  # same, the explicit-multiply path\n",
    "assert x @ y == 32.0, \"dot product 1*4+2*5+3*6 should be 32\"\n",
    "print(\"[ ok ] dot product = 32\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0200f18",
   "metadata": {},
   "source": [
    "> **Interpretation.** Three syntaxes, one number. The `@` operator is the one to reach for; `(x*y).sum()` is the same computation written out, which is occasionally useful when you want to see the per-element contributions before summing.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a191d3c",
   "metadata": {},
   "source": [
    "### Matrix multiplication and the shape rule\n",
    "\n",
    "For $A \\in \\mathbb{R}^{m \\times k}$ and $B \\in \\mathbb{R}^{k \\times n}$, the product $AB \\in \\mathbb{R}^{m \\times n}$ has entries $(AB)_{ij} = \\sum_p A_{ip} B_{pj}$. The rule to burn in: **inner dimensions must match; outer dimensions become the output shape.** `(3, 4) @ (4, 5) = (3, 5)`. `(3, 4) @ (5, 4)` is an error. We trace it on a tiny hand-checkable example first.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "b7a28a3a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.339443Z",
     "iopub.status.busy": "2026-06-10T18:44:02.339381Z",
     "iopub.status.idle": "2026-06-10T18:44:02.341464Z",
     "shell.execute_reply": "2026-06-10T18:44:02.341127Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[19. 22.]\n",
      " [43. 50.]]\n",
      "[ ok ] 2x2 matmul matches hand computation\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1., 2.],\n",
    "              [3., 4.]])\n",
    "Bm = np.array([[5., 6.],\n",
    "               [7., 8.]])\n",
    "# By hand: (AB)[0,0] = 1*5 + 2*7 = 19 ; (AB)[0,1] = 1*6 + 2*8 = 22\n",
    "#          (AB)[1,0] = 3*5 + 4*7 = 43 ; (AB)[1,1] = 3*6 + 4*8 = 50\n",
    "print(A @ Bm)\n",
    "assert np.array_equal(A @ Bm, [[19., 22.], [43., 50.]]), \\\n",
    "    \"2x2 matmul should give [[19,22],[43,50]] by the row-dot-column rule\"\n",
    "print(\"[ ok ] 2x2 matmul matches hand computation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83100d46",
   "metadata": {},
   "source": [
    "> **Common confusion:** the inner dimensions are the ones that touch in the middle. In `(3, 4) @ (4, 5)`, the two 4s touch and vanish; the surviving 3 and 5 are the output shape. If the touching numbers differ, NumPy raises a shape error rather than guessing.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd9a745d",
   "metadata": {},
   "source": [
    "### Exercise 0.2 — Matrix multiply from scratch\n",
    "`Difficulty 2/5 · ~12 min`\n",
    "\n",
    "Implement `matmul(A, B)` with explicit loops over the output indices, returning the same thing `A @ B` does. No `@`, no `np.dot`, no `np.matmul` inside. The check first verifies the shape rule rejects a bad pair, then compares your output against NumPy's `@` on random matrices: NumPy is the oracle, your code is right iff it reproduces it.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "47b991dd",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.342264Z",
     "iopub.status.busy": "2026-06-10T18:44:02.342203Z",
     "iopub.status.idle": "2026-06-10T18:44:02.345728Z",
     "shell.execute_reply": "2026-06-10T18:44:02.345443Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ ok ] 0.2 shape rule\n",
      "[ -- ] 0.2 vs numpy: not attempted yet — fill in the TODO above, then re-run.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def matmul(A, B):\n",
    "    \"\"\"Multiply A (m, k) by B (k, n) -> (m, n) with explicit loops.\"\"\"\n",
    "    A = np.asarray(A, dtype=float)\n",
    "    B = np.asarray(B, dtype=float)\n",
    "    m, k = A.shape\n",
    "    k2, n = B.shape\n",
    "    assert k == k2, f\"inner dims must match: A is (..,{k}), B is ({k2},..)\"\n",
    "    out = np.zeros((m, n))\n",
    "    # TODO 1: triple loop. out[i, j] = sum over p of A[i, p] * B[p, j].\n",
    "    #         (one pedagogical exception to \"never loop\": you are showing what @ does.)\n",
    "    raise NotImplementedError\n",
    "\n",
    "def _matmul_rejects_bad_shapes():\n",
    "    try:\n",
    "        matmul(np.ones((2, 3)), np.ones((4, 5)))\n",
    "    except AssertionError:\n",
    "        return\n",
    "    raise AssertionError(\"matmul should reject (2,3) @ (4,5): inner dims 3 and 4 disagree\")\n",
    "\n",
    "def _matmul_matches_numpy():\n",
    "    P = rng.standard_normal((4, 6)); Q = rng.standard_normal((6, 3))\n",
    "    check_close(matmul(P, Q), P @ Q, msg=\"your matmul disagrees with NumPy's @\")\n",
    "\n",
    "check(\"0.2 shape rule\", _matmul_rejects_bad_shapes)\n",
    "check(\"0.2 vs numpy\",   _matmul_matches_numpy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f7be910",
   "metadata": {},
   "source": [
    "<details><summary>Hint 1 (conceptual)</summary>`out[i, j]` is the dot product of row `i` of `A` with column `j` of `B`. Three nested loops: `i` over `m`, `j` over `n`, `p` over `k`.</details>\n",
    "\n",
    "<details><summary>Hint 2 (pseudocode)</summary>\n",
    "\n",
    "```python\n",
    "for i in range(m):\n",
    "    for j in range(n):\n",
    "        s = 0.0\n",
    "        for p in range(k):\n",
    "            s += A[i, p] * B[p, j]\n",
    "        out[i, j] = s\n",
    "return out\n",
    "```\n",
    "</details>\n",
    "\n",
    "<details><summary>Help — \"shapes wrong\" or values off by transpose</summary>If your result is `(n, m)` instead of `(m, n)`, you swapped the index order: the row index `i` selects from `A` and indexes the output's first axis. If values are wrong but the shape is right, check that the contracted index `p` ranges over `k` and indexes `A`'s second axis and `B`'s first.</details>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "4c8aac84",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.346392Z",
     "iopub.status.busy": "2026-06-10T18:44:02.346325Z",
     "iopub.status.idle": "2026-06-10T18:44:02.349632Z",
     "shell.execute_reply": "2026-06-10T18:44:02.349369Z"
    },
    "collapsed": true,
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ ok ] 0.2 shape rule\n",
      "[ ok ] 0.2 vs numpy\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#@title Solution { display-mode: \"form\" }\n",
    "# solution: redefines matmul; the checks below re-verify the reference.\n",
    "def matmul(A, B):\n",
    "    A = np.asarray(A, dtype=float)\n",
    "    B = np.asarray(B, dtype=float)\n",
    "    m, k = A.shape\n",
    "    k2, n = B.shape\n",
    "    assert k == k2, f\"inner dims must match: A is (..,{k}), B is ({k2},..)\"\n",
    "    out = np.zeros((m, n))\n",
    "    for i in range(m):\n",
    "        for j in range(n):\n",
    "            s = 0.0\n",
    "            for p in range(k):\n",
    "                s += A[i, p] * B[p, j]\n",
    "            out[i, j] = s\n",
    "    return out\n",
    "\n",
    "check(\"0.2 shape rule\", _matmul_rejects_bad_shapes, required=True)\n",
    "check(\"0.2 vs numpy\",   _matmul_matches_numpy, required=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13e7738c",
   "metadata": {},
   "source": [
    "> **Interpretation.** The triple loop is exactly what `@` does, just slower by a few hundred times because NumPy dispatches to a compiled BLAS kernel. You will never write this loop in real code. Writing it once is how you stop being scared of the operation.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d94852a8",
   "metadata": {},
   "source": [
    "### Norms measure size\n",
    "\n",
    "The $L^2$ (Euclidean) norm $\\|\\mathbf{x}\\|_2 = \\sqrt{\\sum_i x_i^2}$ is a length; you use it for distances and energies. The $L^1$ norm $\\|\\mathbf{x}\\|_1 = \\sum_i |x_i|$ shows up as a sparsity penalty. The Frobenius norm of a matrix is the $L^2$ norm of its flattened entries. NumPy spells all of them `np.linalg.norm`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "60e410e1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.350429Z",
     "iopub.status.busy": "2026-06-10T18:44:02.350366Z",
     "iopub.status.idle": "2026-06-10T18:44:02.352554Z",
     "shell.execute_reply": "2026-06-10T18:44:02.352170Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L2: 5.0\n",
      "L1: 7.0\n",
      "Linf: 4.0\n",
      "[ ok ] norms of [3,4]\n"
     ]
    }
   ],
   "source": [
    "v = np.array([3., 4.])\n",
    "print(\"L2:\", np.linalg.norm(v))           # sqrt(9+16) = 5, the 3-4-5 triangle\n",
    "print(\"L1:\", np.linalg.norm(v, ord=1))    # |3| + |4| = 7\n",
    "print(\"Linf:\", np.linalg.norm(v, ord=np.inf))  # max(|3|,|4|) = 4\n",
    "assert np.isclose(np.linalg.norm(v), 5.0), \"L2 norm of [3,4] is the 3-4-5 hypotenuse\"\n",
    "print(\"[ ok ] norms of [3,4]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecac8721",
   "metadata": {},
   "source": [
    "> **Key takeaways.** The dot product is the workhorse scalar. Matmul is row-dot-column with the inner-dimension shape rule. Norms measure size; reach for `np.linalg.norm` with an `ord` argument rather than rewriting the formula. Your hand-written matmul reproduced `@` to floating-point tolerance, which is the proof you understand it.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee74939a",
   "metadata": {},
   "source": [
    "## Part 3 — Calculus and the gradient\n",
    "\n",
    "> **Objectives.** Estimate a derivative numerically and watch it converge to the analytic one; apply the chain rule; derive the least-squares gradient and confirm it against `torch.autograd`.\n",
    "\n",
    "We open the way Karpathy opens micrograd: with the definition of a derivative, made into running code. The derivative is the limit of a slope,\n",
    "$$\\frac{df}{dx} = \\lim_{h \\to 0} \\frac{f(x + h) - f(x)}{h}.$$\n",
    "We cannot take a limit on a computer, but we can take a small `h` and see what the slope approaches.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "6a160816",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.353575Z",
     "iopub.status.busy": "2026-06-10T18:44:02.353466Z",
     "iopub.status.idle": "2026-06-10T18:44:02.356585Z",
     "shell.execute_reply": "2026-06-10T18:44:02.356126Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "h=1e-01  forward-difference slope = 5.611000\n",
      "h=1e-03  forward-difference slope = 4.510506\n",
      "h=1e-05  forward-difference slope = 4.500105\n",
      "h=1e-07  forward-difference slope = 4.500001\n",
      "analytic 4x^3 - 6x at x=1.5: 4.500000\n"
     ]
    }
   ],
   "source": [
    "def f(x):\n",
    "    return x**4 - 3*x**2 + 2   # analytic derivative: 4*x**3 - 6*x\n",
    "\n",
    "x0 = 1.5\n",
    "for h in [1e-1, 1e-3, 1e-5, 1e-7]:\n",
    "    numerical = (f(x0 + h) - f(x0)) / h\n",
    "    print(f\"h={h:.0e}  forward-difference slope = {numerical:.6f}\")\n",
    "analytic = 4*x0**3 - 6*x0\n",
    "print(f\"analytic 4x^3 - 6x at x={x0}: {analytic:.6f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c80a999",
   "metadata": {},
   "source": [
    "> **Interpretation.** As `h` shrinks the forward-difference estimate marches toward the analytic value `4*1.5**3 - 6*1.5 = 4.5`. It never lands exactly: too large an `h` and the slope is curved over the interval; too small and floating-point subtraction loses precision. This tension is the whole reason `h=1e-5` is a common default rather than `1e-12`.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcfd53ca",
   "metadata": {},
   "source": [
    "### A centered difference is more accurate\n",
    "\n",
    "The forward difference uses `f(x+h) - f(x)`. The centered difference uses points on both sides, `(f(x+h) - f(x-h)) / (2h)`, which cancels the leading error term and is accurate to order $h^2$ instead of $h$. This is the finite-difference gradient check you will use to validate hand-derived gradients for the rest of the curriculum.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "3140b4db",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.357367Z",
     "iopub.status.busy": "2026-06-10T18:44:02.357295Z",
     "iopub.status.idle": "2026-06-10T18:44:02.359484Z",
     "shell.execute_reply": "2026-06-10T18:44:02.359012Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "centered estimate: 4.50000000067341\n",
      "analytic:          4.5\n",
      "[ ok ] centered difference matches analytic derivative\n"
     ]
    }
   ],
   "source": [
    "def numerical_derivative(fn, x, h=1e-5):\n",
    "    return (fn(x + h) - fn(x - h)) / (2 * h)   # centered difference\n",
    "\n",
    "print(\"centered estimate:\", numerical_derivative(f, x0))\n",
    "print(\"analytic:         \", analytic)\n",
    "assert abs(numerical_derivative(f, x0) - analytic) < 1e-6, \\\n",
    "    \"centered difference should match the analytic derivative to ~1e-6\"\n",
    "print(\"[ ok ] centered difference matches analytic derivative\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cef500d0",
   "metadata": {},
   "source": [
    "### Exercise 0.3 — A finite-difference gradient checker\n",
    "`Difficulty 3/5 · ~15 min`\n",
    "\n",
    "This is the workhorse you will reuse whenever you derive a gradient by hand. Fill in `grad_check(fn, grad_fn, x)` so it returns the absolute difference between the analytic gradient `grad_fn(x)` and a centered finite-difference estimate of `fn`'s derivative at `x`. The check feeds it a function whose gradient you have on paper, asserts the gap is tiny, and then asserts that a deliberately WRONG analytic gradient is caught (a checker that never complains is useless).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f2a08030",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.360263Z",
     "iopub.status.busy": "2026-06-10T18:44:02.360199Z",
     "iopub.status.idle": "2026-06-10T18:44:02.363936Z",
     "shell.execute_reply": "2026-06-10T18:44:02.363735Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ -- ] 0.3 catches correct: not attempted yet — fill in the TODO above, then re-run.\n",
      "[ -- ] 0.3 catches wrong: not attempted yet — fill in the TODO above, then re-run.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def grad_check(fn, grad_fn, x, h=1e-5):\n",
    "    \"\"\"|analytic grad - centered finite difference| at scalar x.\"\"\"\n",
    "    # TODO 1: centered finite difference of fn at x with step h\n",
    "    numerical = None\n",
    "    # TODO 2: the analytic gradient the caller supplied\n",
    "    analytic = None\n",
    "    attempted(numerical, analytic)\n",
    "    # TODO 3: return the absolute difference\n",
    "    return abs(numerical - analytic)\n",
    "\n",
    "def _grad_check_passes_correct():\n",
    "    g = lambda x: x**3            # f\n",
    "    dg = lambda x: 3 * x**2       # correct f'\n",
    "    assert grad_check(g, dg, 2.0) < 1e-4, \\\n",
    "        \"with the correct gradient the gap should be ~1e-6, not large\"\n",
    "\n",
    "def _grad_check_catches_wrong():\n",
    "    g = lambda x: x**3\n",
    "    wrong = lambda x: 2 * x**2    # WRONG: off by a factor and a power\n",
    "    assert grad_check(g, wrong, 2.0) > 1.0, \\\n",
    "        \"a wrong analytic gradient must produce a large gap; the checker would be useless otherwise\"\n",
    "\n",
    "check(\"0.3 catches correct\", _grad_check_passes_correct)\n",
    "check(\"0.3 catches wrong\",   _grad_check_catches_wrong)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d5919a2",
   "metadata": {},
   "source": [
    "<details><summary>Hint 1 (conceptual)</summary>You already wrote the centered difference in `numerical_derivative` above. The analytic value is whatever `grad_fn(x)` returns. The result is `abs` of their difference.</details>\n",
    "\n",
    "<details><summary>Hint 2 (pseudocode)</summary>\n",
    "\n",
    "```python\n",
    "numerical = (fn(x + h) - fn(x - h)) / (2 * h)\n",
    "analytic  = grad_fn(x)\n",
    "return abs(numerical - analytic)\n",
    "```\n",
    "</details>\n",
    "\n",
    "<details><summary>Help — \"the wrong-gradient check fails\"</summary>If `_grad_check_catches_wrong` fails, your function is probably returning the numerical value alone and ignoring `grad_fn`, so it never disagrees with anything. Make sure you actually subtract `grad_fn(x)`.</details>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "33fa1454",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.364678Z",
     "iopub.status.busy": "2026-06-10T18:44:02.364620Z",
     "iopub.status.idle": "2026-06-10T18:44:02.367068Z",
     "shell.execute_reply": "2026-06-10T18:44:02.366735Z"
    },
    "collapsed": true,
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ ok ] 0.3 catches correct\n",
      "[ ok ] 0.3 catches wrong\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#@title Solution { display-mode: \"form\" }\n",
    "# solution: redefines grad_check; the checks below re-verify the reference.\n",
    "def grad_check(fn, grad_fn, x, h=1e-5):\n",
    "    numerical = (fn(x + h) - fn(x - h)) / (2 * h)\n",
    "    analytic = grad_fn(x)\n",
    "    return abs(numerical - analytic)\n",
    "\n",
    "check(\"0.3 catches correct\", _grad_check_passes_correct, required=True)\n",
    "check(\"0.3 catches wrong\",   _grad_check_catches_wrong, required=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1427ac15",
   "metadata": {},
   "source": [
    "> **Interpretation.** A gradient checker is only as good as its ability to fail. The second assert, that a wrong gradient produces a gap above 1.0, is the part most people forget. Without it you have a checker that rubber-stamps anything.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc69f366",
   "metadata": {},
   "source": [
    "### The chain rule\n",
    "\n",
    "If $y = f(u)$ and $u = g(x)$, then $\\frac{dy}{dx} = \\frac{dy}{du}\\cdot\\frac{du}{dx}$. This is the rule backpropagation is built on. Take $y = (3x + 2)^4$: with $u = 3x + 2$, we get $\\frac{dy}{du} = 4u^3$ and $\\frac{du}{dx} = 3$, so $\\frac{dy}{dx} = 12(3x + 2)^3$. We confirm the hand answer with the checker we just built.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "6a935582",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.367871Z",
     "iopub.status.busy": "2026-06-10T18:44:02.367801Z",
     "iopub.status.idle": "2026-06-10T18:44:02.369931Z",
     "shell.execute_reply": "2026-06-10T18:44:02.369612Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gap between hand chain-rule gradient and finite difference: 4.80e-08\n",
      "[ ok ] chain rule on (3x+2)^4\n"
     ]
    }
   ],
   "source": [
    "composite = lambda x: (3*x + 2)**4\n",
    "chain_grad = lambda x: 12 * (3*x + 2)**3   # dy/du * du/dx = 4u^3 * 3\n",
    "gap = grad_check(composite, chain_grad, 0.7)\n",
    "print(f\"gap between hand chain-rule gradient and finite difference: {gap:.2e}\")\n",
    "assert gap < 1e-5, \"chain-rule gradient 12(3x+2)^3 should match the numerical derivative\"\n",
    "print(\"[ ok ] chain rule on (3x+2)^4\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8ef2fc2",
   "metadata": {},
   "source": [
    "### The gradient that runs the whole curriculum\n",
    "\n",
    "For the least-squares loss $L = \\tfrac{1}{2}\\|X\\mathbf{w} - \\mathbf{y}\\|^2$, the gradient is\n",
    "$$\\nabla_{\\mathbf{w}} L = X^\\top (X\\mathbf{w} - \\mathbf{y}).$$\n",
    "Every variable in that formula has a name in the code below. When in doubt, never trust a hand-derived gradient until `torch.autograd` agrees with it. This vector-valued check is the multi-dimensional version of the scalar `grad_check` above.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "8ff8d2ba",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.370663Z",
     "iopub.status.busy": "2026-06-10T18:44:02.370595Z",
     "iopub.status.idle": "2026-06-10T18:44:02.617213Z",
     "shell.execute_reply": "2026-06-10T18:44:02.616697Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "autograd w.grad: [ 0.3108 12.6724 21.067 ]\n",
      "hand-derived:    [ 0.3108 12.6724 21.067 ]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ ok ] hand gradient X^T(Xw - y) matches torch.autograd\n"
     ]
    }
   ],
   "source": [
    "torch.manual_seed(SEED)\n",
    "Xt = torch.randn(10, 3)\n",
    "yt = torch.randn(10)\n",
    "w = torch.randn(3, requires_grad=True)\n",
    "\n",
    "loss = ((Xt @ w - yt) ** 2).sum() / 2   # L = 1/2 ||Xw - y||^2\n",
    "loss.backward()                          # autograd fills w.grad\n",
    "\n",
    "with torch.no_grad():\n",
    "    manual = Xt.T @ (Xt @ w - yt)        # hand-derived: X^T (Xw - y)\n",
    "print(\"autograd w.grad:\", w.grad.numpy().round(4))\n",
    "print(\"hand-derived:   \", manual.numpy().round(4))\n",
    "torch.testing.assert_close(w.grad, manual, atol=1e-5, rtol=1e-4)\n",
    "print(\"[ ok ] hand gradient X^T(Xw - y) matches torch.autograd\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e49dc00",
   "metadata": {},
   "source": [
    "> **Note:** we wrap the hand computation in `torch.no_grad()` because it is a check, not part of any graph we will differentiate; that is the inference-hygiene habit the deep-learning chapters lean on. If the `assert_close` had failed, the math would be wrong, not the library.\n",
    "\n",
    "> **Key takeaways.** A derivative is a limit of a slope; on a computer you approximate it with a small step, centered for accuracy. The chain rule multiplies local derivatives along the composition, which is all backprop is. Always grad-check: a hand gradient is a hypothesis until autograd confirms it.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b17b96df",
   "metadata": {},
   "source": [
    "## Part 4 — Probability and maximum likelihood\n",
    "\n",
    "> **Objectives.** Compute expectation and variance from a distribution; apply Bayes' rule on a concrete table; understand why negative log-likelihood is the loss behind nearly every model.\n",
    "\n",
    "The expected value of a discrete random variable is $E[X] = \\sum_x x\\, p(x)$, the probability-weighted average. The variance is $\\operatorname{Var}(X) = E[X^2] - E[X]^2$. We sanity-check the formula against the law of large numbers: draw many samples and watch the empirical mean approach the analytic one.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "2fec0568",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.618292Z",
     "iopub.status.busy": "2026-06-10T18:44:02.618156Z",
     "iopub.status.idle": "2026-06-10T18:44:02.621955Z",
     "shell.execute_reply": "2026-06-10T18:44:02.621551Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "analytic  E[X]=3.5000  Var(X)=2.9167\n",
      "empirical E[X]=3.4980  Var(X)=2.9103\n",
      "[ ok ] empirical mean approaches the analytic expectation\n"
     ]
    }
   ],
   "source": [
    "rng = np.random.default_rng(SEED)   # re-seed so this cell reproduces in isolation\n",
    "values = np.array([1, 2, 3, 4, 5, 6])           # a fair die\n",
    "probs = np.full(6, 1/6)\n",
    "EX = (values * probs).sum()                      # 3.5\n",
    "VarX = (values**2 * probs).sum() - EX**2         # E[X^2] - E[X]^2\n",
    "print(f\"analytic  E[X]={EX:.4f}  Var(X)={VarX:.4f}\")\n",
    "\n",
    "samples = rng.integers(1, 7, size=100_000)       # 100k die rolls\n",
    "print(f\"empirical E[X]={samples.mean():.4f}  Var(X)={samples.var():.4f}\")\n",
    "assert abs(samples.mean() - EX) < 0.05, \"empirical mean of 100k rolls should be near 3.5\"\n",
    "print(\"[ ok ] empirical mean approaches the analytic expectation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aed0fade",
   "metadata": {},
   "source": [
    "> **Interpretation.** The analytic $E[X] = 3.5$ and the average of 100,000 rolls agree to two decimals. This is the law of large numbers doing what we will lean on constantly: a sample mean is an estimate of an expectation, and it gets better with more samples. Note we used `E[X^2] - E[X]^2` for variance here because the distribution is exact; on near-equal floating-point scales that formula can cancel catastrophically, which the safety lens revisits.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ed70544",
   "metadata": {},
   "source": [
    "### Exercise 0.4 — Bayes' rule on a medical test\n",
    "`Difficulty 2/5 · ~10 min`\n",
    "\n",
    "A disease has prevalence $p(D) = 0.01$. A test has sensitivity $p(+\\mid D) = 0.99$ and false-positive rate $p(+\\mid \\neg D) = 0.05$. Fill in `posterior(...)` to return $p(D \\mid +)$ via Bayes' rule. The check asserts the famous counterintuitive answer (it is far below 0.99) and that the posterior of a perfectly specific test ($p(+\\mid\\neg D)=0$) is exactly 1.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "001f03d1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.623065Z",
     "iopub.status.busy": "2026-06-10T18:44:02.622990Z",
     "iopub.status.idle": "2026-06-10T18:44:02.627239Z",
     "shell.execute_reply": "2026-06-10T18:44:02.626888Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ -- ] 0.4 rare-disease posterior: not attempted yet — fill in the TODO above, then re-run.\n",
      "[ -- ] 0.4 perfect specificity: not attempted yet — fill in the TODO above, then re-run.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def posterior(p_d, p_pos_given_d, p_pos_given_notd):\n",
    "    \"\"\"p(D | +) by Bayes' rule. All inputs in [0, 1].\"\"\"\n",
    "    # TODO 1: evidence p(+) = p(+|D)p(D) + p(+|~D)p(~D)\n",
    "    evidence = None\n",
    "    attempted(evidence)\n",
    "    # TODO 2: posterior = p(+|D) p(D) / p(+)\n",
    "    return (p_pos_given_d * p_d) / evidence\n",
    "\n",
    "def _bayes_counterintuitive():\n",
    "    post = posterior(0.01, 0.99, 0.05)\n",
    "    check_close(post, 0.99 * 0.01 / (0.99 * 0.01 + 0.05 * 0.99), atol=1e-9,\n",
    "                msg=\"p(D|+) for a rare disease is far below the 0.99 sensitivity\")\n",
    "    assert post < 0.2, f\"with 1% prevalence the posterior should be ~0.17, got {post:.3f}\"\n",
    "\n",
    "def _bayes_perfect_specificity():\n",
    "    assert abs(posterior(0.01, 0.99, 0.0) - 1.0) < 1e-12, \\\n",
    "        \"if there are no false positives, a positive test is certain disease\"\n",
    "\n",
    "check(\"0.4 rare-disease posterior\", _bayes_counterintuitive)\n",
    "check(\"0.4 perfect specificity\",    _bayes_perfect_specificity)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ba75778",
   "metadata": {},
   "source": [
    "<details><summary>Hint 1 (conceptual)</summary>The denominator is the total probability of a positive test, summing the two ways it can happen: a true positive and a false positive. $p(\\neg D) = 1 - p(D)$.</details>\n",
    "\n",
    "<details><summary>Hint 2 (pseudocode)</summary>\n",
    "\n",
    "```python\n",
    "evidence = p_pos_given_d * p_d + p_pos_given_notd * (1 - p_d)\n",
    "return p_pos_given_d * p_d / evidence\n",
    "```\n",
    "</details>\n",
    "\n",
    "<details><summary>Help — \"posterior is way too high, near 0.99\"</summary>You likely forgot the false-positive term in the evidence, or used `p_d` where you needed `1 - p_d`. With 99 healthy people per 1 sick person, the 5% false-positive rate produces far more positives than the disease does.</details>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "31d7d461",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.628068Z",
     "iopub.status.busy": "2026-06-10T18:44:02.627994Z",
     "iopub.status.idle": "2026-06-10T18:44:02.630598Z",
     "shell.execute_reply": "2026-06-10T18:44:02.630261Z"
    },
    "collapsed": true,
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ ok ] 0.4 rare-disease posterior\n",
      "[ ok ] 0.4 perfect specificity\n",
      "p(D | +) = 0.167\n"
     ]
    }
   ],
   "source": [
    "#@title Solution { display-mode: \"form\" }\n",
    "# solution: redefines posterior; the checks below re-verify the reference.\n",
    "def posterior(p_d, p_pos_given_d, p_pos_given_notd):\n",
    "    evidence = p_pos_given_d * p_d + p_pos_given_notd * (1 - p_d)\n",
    "    return (p_pos_given_d * p_d) / evidence\n",
    "\n",
    "check(\"0.4 rare-disease posterior\", _bayes_counterintuitive, required=True)\n",
    "check(\"0.4 perfect specificity\",    _bayes_perfect_specificity, required=True)\n",
    "print(f\"p(D | +) = {posterior(0.01, 0.99, 0.05):.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecae5aef",
   "metadata": {},
   "source": [
    "> **Interpretation.** Even with a 99%-sensitive test, a positive result means only about a 17% chance of disease, because the disease is rare and false positives outnumber true ones. This is the base-rate fallacy, and it is the same arithmetic as a \"99% accurate\" classifier on a rare-positive stream, which Ch 03 takes apart.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5478317",
   "metadata": {},
   "source": [
    "### Maximum likelihood is why losses look the way they do\n",
    "\n",
    "Given data and a model $p_\\theta(y\\mid x)$, the maximum-likelihood estimate maximizes $\\prod_i p_\\theta(y_i\\mid x_i)$. Products of probabilities underflow, so we maximize the log instead, equivalently minimizing the **negative log-likelihood** $-\\sum_i \\log p_\\theta(y_i\\mid x_i)$. For a Gaussian model this is mean squared error; for a categorical model it is cross-entropy. We show the underflow that forces the log, on a long sequence of small probabilities.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "9814a7d6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.631424Z",
     "iopub.status.busy": "2026-06-10T18:44:02.631347Z",
     "iopub.status.idle": "2026-06-10T18:44:02.633559Z",
     "shell.execute_reply": "2026-06-10T18:44:02.633262Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "raw product of 2000 probs: 0.0   (underflowed to 0)\n",
      "sum of log-probs:          -2407.95   (a usable number)\n",
      "[ ok ] products underflow, log-sums do not\n"
     ]
    }
   ],
   "source": [
    "ps = np.full(2000, 0.3)              # 2000 independent events each with p=0.3\n",
    "prod = np.prod(ps)                   # the raw likelihood\n",
    "log_sum = np.log(ps).sum()           # the log-likelihood\n",
    "print(f\"raw product of 2000 probs: {prod}   (underflowed to 0)\")\n",
    "print(f\"sum of log-probs:          {log_sum:.2f}   (a usable number)\")\n",
    "assert prod == 0.0 and log_sum < 0, \"the raw product underflows to 0; the log-sum stays finite\"\n",
    "print(\"[ ok ] products underflow, log-sums do not\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0dd44cde",
   "metadata": {},
   "source": [
    "> **Key takeaways.** Expectation is a probability-weighted average; a sample mean estimates it. Bayes' rule flips a conditional and the base rate dominates rare events. Maximum likelihood minimizes negative log-likelihood; we take logs because products of many probabilities underflow to exactly zero. MSE and cross-entropy are both NLL under different noise assumptions.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a01c062d",
   "metadata": {},
   "source": [
    "## Part 5 — Information theory\n",
    "\n",
    "> **Objectives.** Compute entropy and cross-entropy; build `softmax` and `cross_entropy` from scratch and check them against PyTorch; experience the naive-softmax overflow bug and fix it.\n",
    "\n",
    "Entropy measures uncertainty: $H(p) = -\\sum_x p(x)\\log p(x)$, in nats when the log is natural. A certain outcome has entropy 0; a uniform distribution over $K$ outcomes has the maximum, $\\log K$. We compute the fair-coin and certain-coin entropies from the prediction at the top of the notebook.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "0a7d742d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.634444Z",
     "iopub.status.busy": "2026-06-10T18:44:02.634369Z",
     "iopub.status.idle": "2026-06-10T18:44:02.637105Z",
     "shell.execute_reply": "2026-06-10T18:44:02.636655Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fair coin  H = 0.6931 nats  (= ln 2 = 0.6931)\n",
      "certain    H = -0.0000 nats\n",
      "uniform-4  H = 1.3863 nats  (= ln 4 = 1.3863)\n",
      "[ ok ] entropy of fair / certain / uniform\n"
     ]
    }
   ],
   "source": [
    "def entropy(p):\n",
    "    p = np.asarray(p, dtype=float)\n",
    "    p = p[p > 0]                       # 0*log0 := 0; drop zeros to avoid log(0)\n",
    "    return float(-(p * np.log(p)).sum())\n",
    "\n",
    "print(f\"fair coin  H = {entropy([0.5, 0.5]):.4f} nats  (= ln 2 = {np.log(2):.4f})\")\n",
    "print(f\"certain    H = {entropy([1.0, 0.0]):.4f} nats\")\n",
    "print(f\"uniform-4  H = {entropy([0.25]*4):.4f} nats  (= ln 4 = {np.log(4):.4f})\")\n",
    "assert np.isclose(entropy([0.5, 0.5]), np.log(2)), \"fair-coin entropy is ln 2\"\n",
    "assert entropy([1.0, 0.0]) == 0.0, \"a certain outcome carries no surprise\"\n",
    "print(\"[ ok ] entropy of fair / certain / uniform\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5eff27cc",
   "metadata": {},
   "source": [
    "> **Interpretation.** The fair coin sits at `ln 2 ≈ 0.693` nats, the certain coin at 0, and the 4-way uniform at `ln 4`. Entropy is the average surprise of a sample; a distribution you can predict perfectly has none.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3aa4da9",
   "metadata": {},
   "source": [
    "### Cross-entropy is the classifier loss\n",
    "\n",
    "Cross-entropy between truth $p$ and model $q$ is $H(p, q) = -\\sum_x p(x)\\log q(x)$. When $p$ is a one-hot target, this collapses to $-\\log q(\\text{true class})$, which is exactly the negative log-likelihood. PyTorch's `F.cross_entropy` takes raw logits and a class index and does `log_softmax` then `nll_loss` internally. We confirm the two paths agree.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "48d836e1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.638030Z",
     "iopub.status.busy": "2026-06-10T18:44:02.637961Z",
     "iopub.status.idle": "2026-06-10T18:44:02.641756Z",
     "shell.execute_reply": "2026-06-10T18:44:02.641436Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F.cross_entropy: 1.130354\n",
      "manual NLL:      1.130354\n",
      "[ ok ] F.cross_entropy == -mean log-softmax at the target = NLL\n"
     ]
    }
   ],
   "source": [
    "import torch.nn.functional as F\n",
    "torch.manual_seed(SEED)\n",
    "logits = torch.randn(4, 3)               # (batch, n_classes)\n",
    "targets = torch.tensor([0, 2, 1, 1])     # integer class labels\n",
    "ce = F.cross_entropy(logits, targets)\n",
    "\n",
    "log_probs = F.log_softmax(logits, dim=-1)             # log of the softmax\n",
    "manual = -log_probs[torch.arange(4), targets].mean()  # pick the true-class log-prob, negate, average\n",
    "print(f\"F.cross_entropy: {ce.item():.6f}\")\n",
    "print(f\"manual NLL:      {manual.item():.6f}\")\n",
    "torch.testing.assert_close(ce, manual, atol=1e-6, rtol=1e-5)\n",
    "print(\"[ ok ] F.cross_entropy == -mean log-softmax at the target = NLL\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7cb3dea",
   "metadata": {},
   "source": [
    "> **Interpretation.** Cross-entropy and negative log-likelihood are two names for the same number when the target is a one-hot. `F.cross_entropy` is the numerically stable, fused version; the manual `log_softmax` then gather then negate then mean is what it computes.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de497cf3",
   "metadata": {},
   "source": [
    "### A deliberate failure: the softmax that returns `nan`\n",
    "\n",
    "Softmax is $q_k = e^{z_k} / \\sum_j e^{z_j}$. Translated to code literally, it overflows: `exp` of a large logit is `inf`, and `inf / inf` is `nan`. We run the naive version on a large logit on purpose, watch it fail, then fix it with the log-sum-exp trick (subtract the max before exponentiating, which is exact because it cancels in the ratio). Staging the bug is the point; this exact footgun has cost real training runs.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "dfab1c0b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.642663Z",
     "iopub.status.busy": "2026-06-10T18:44:02.642589Z",
     "iopub.status.idle": "2026-06-10T18:44:02.645112Z",
     "shell.execute_reply": "2026-06-10T18:44:02.644845Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "naive softmax on a large logit: [ 0.  0. nan]\n",
      "contains nan? True\n",
      "[ ok ] reproduced the overflow bug on purpose\n"
     ]
    }
   ],
   "source": [
    "def softmax_naive(z):\n",
    "    z = np.asarray(z, dtype=float)\n",
    "    e = np.exp(z)                 # overflows for large z\n",
    "    return e / e.sum()\n",
    "\n",
    "bad = softmax_naive([1.0, 2.0, 1000.0])   # 1000 is plausible for an unnormalized logit\n",
    "print(\"naive softmax on a large logit:\", bad)\n",
    "print(\"contains nan?\", np.isnan(bad).any())\n",
    "assert np.isnan(bad).any(), \"the naive softmax SHOULD produce nan here; that is the bug we are demonstrating\"\n",
    "print(\"[ ok ] reproduced the overflow bug on purpose\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5fe37a7f",
   "metadata": {},
   "source": [
    "> **Common confusion:** the fix is not to clip or to use higher precision. Subtracting the maximum logit before exponentiating leaves the ratio unchanged (the constant cancels top and bottom) but caps the largest exponent at `exp(0) = 1`, so nothing overflows. This is the log-sum-exp trick, and it is why you should call `F.cross_entropy` rather than write softmax-then-log by hand.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "8c80df54",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.645952Z",
     "iopub.status.busy": "2026-06-10T18:44:02.645881Z",
     "iopub.status.idle": "2026-06-10T18:44:02.648784Z",
     "shell.execute_reply": "2026-06-10T18:44:02.648361Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "stable softmax on the same input: [0. 0. 1.]\n",
      "[ ok ] stable softmax: no nan, sums to 1, argmax preserved\n"
     ]
    }
   ],
   "source": [
    "def softmax_stable(z):\n",
    "    z = np.asarray(z, dtype=float)\n",
    "    z = z - z.max()               # subtract max: ratio unchanged, largest exp is 1\n",
    "    e = np.exp(z)\n",
    "    return e / e.sum()\n",
    "\n",
    "good = softmax_stable([1.0, 2.0, 1000.0])\n",
    "print(\"stable softmax on the same input:\", good)\n",
    "assert not np.isnan(good).any(), \"the stable softmax must not produce nan\"\n",
    "assert np.isclose(good.sum(), 1.0), \"softmax outputs are a distribution and sum to 1\"\n",
    "assert good.argmax() == 2, \"the largest logit should still get the most mass\"\n",
    "print(\"[ ok ] stable softmax: no nan, sums to 1, argmax preserved\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68b259b5",
   "metadata": {},
   "source": [
    "### Exercise 0.5 — `softmax` and `cross_entropy` from scratch\n",
    "`Difficulty 3/5 · ~18 min`\n",
    "\n",
    "Implement the stable softmax over the last axis of a `(batch, n_classes)` logit array, then `cross_entropy(logits, targets)` returning the mean NLL. The check verifies rows sum to 1, that your softmax survives a huge logit (the bug above), and that your cross-entropy matches `torch.nn.functional.cross_entropy` to tolerance. PyTorch is the oracle.\n",
    "\n",
    "Harder: make `softmax` work for an arbitrary axis via a keyword argument.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "0a48a86c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.649746Z",
     "iopub.status.busy": "2026-06-10T18:44:02.649582Z",
     "iopub.status.idle": "2026-06-10T18:44:02.654787Z",
     "shell.execute_reply": "2026-06-10T18:44:02.654480Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ -- ] 0.5 softmax distribution: not attempted yet — fill in the TODO above, then re-run.\n",
      "[ -- ] 0.5 softmax no overflow: not attempted yet — fill in the TODO above, then re-run.\n",
      "[ -- ] 0.5 cross_entropy vs torch: not attempted yet — fill in the TODO above, then re-run.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def softmax(logits, axis=-1):\n",
    "    \"\"\"Numerically stable softmax over `axis`. logits: (..., K) -> same shape.\"\"\"\n",
    "    logits = np.asarray(logits, dtype=float)\n",
    "    # TODO 1: subtract the max along `axis` (keepdims=True) to avoid overflow\n",
    "    shifted = None\n",
    "    attempted(shifted)\n",
    "    # TODO 2: exponentiate, then divide by the sum along `axis` (keepdims=True)\n",
    "    e = np.exp(shifted)\n",
    "    return e / e.sum(axis=axis, keepdims=True)\n",
    "\n",
    "def cross_entropy(logits, targets):\n",
    "    \"\"\"Mean NLL. logits: (N, K); targets: (N,) integer class indices.\"\"\"\n",
    "    logits = np.asarray(logits, dtype=float)\n",
    "    targets = np.asarray(targets)\n",
    "    probs = softmax(logits, axis=-1)                 # (N, K)\n",
    "    # TODO 3: gather the probability of the true class for each row, then mean of -log\n",
    "    true_class_probs = None\n",
    "    attempted(true_class_probs)\n",
    "    return float(-np.log(true_class_probs).mean())\n",
    "\n",
    "def _softmax_is_a_distribution():\n",
    "    P = softmax(rng.standard_normal((6, 4)), axis=-1)\n",
    "    check_shape(P, (6, 4))\n",
    "    check_close(P.sum(axis=-1), np.ones(6), msg=\"softmax rows must sum to 1\")\n",
    "\n",
    "def _softmax_survives_overflow():\n",
    "    out = softmax(np.array([1.0, 2.0, 1000.0]))\n",
    "    assert not np.isnan(out).any(), \"softmax must subtract the max so a big logit does not overflow\"\n",
    "\n",
    "def _cross_entropy_matches_torch():\n",
    "    L = rng.standard_normal((5, 3)); t = np.array([0, 2, 1, 1, 0])\n",
    "    want = F.cross_entropy(torch.tensor(L), torch.tensor(t)).item()\n",
    "    check_close(cross_entropy(L, t), want, msg=\"your cross_entropy disagrees with F.cross_entropy\")\n",
    "\n",
    "check(\"0.5 softmax distribution\", _softmax_is_a_distribution)\n",
    "check(\"0.5 softmax no overflow\",  _softmax_survives_overflow)\n",
    "check(\"0.5 cross_entropy vs torch\", _cross_entropy_matches_torch)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2753fe1a",
   "metadata": {},
   "source": [
    "<details><summary>Hint 1 (conceptual)</summary>The stable softmax subtracts the per-row max before `exp`. For cross-entropy, you need the predicted probability of the *correct* class in each row, which is fancy indexing: `probs[np.arange(N), targets]`.</details>\n",
    "\n",
    "<details><summary>Hint 2 (pseudocode)</summary>\n",
    "\n",
    "```python\n",
    "shifted = logits - logits.max(axis=axis, keepdims=True)\n",
    "# cross_entropy:\n",
    "N = logits.shape[0]\n",
    "true_class_probs = probs[np.arange(N), targets]\n",
    "return float(-np.log(true_class_probs).mean())\n",
    "```\n",
    "</details>\n",
    "\n",
    "<details><summary>Help — \"cross_entropy off by a constant or sign\"</summary>If your value is the negative of PyTorch's, you forgot the minus on `log`. If it is the sum rather than the mean, you used `.sum()` where `F.cross_entropy` averages over the batch. If only a huge-logit case fails, your softmax is not subtracting the max.</details>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "959b403c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.655881Z",
     "iopub.status.busy": "2026-06-10T18:44:02.655801Z",
     "iopub.status.idle": "2026-06-10T18:44:02.659418Z",
     "shell.execute_reply": "2026-06-10T18:44:02.659177Z"
    },
    "collapsed": true,
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ ok ] 0.5 softmax distribution\n",
      "[ ok ] 0.5 softmax no overflow\n",
      "[ ok ] 0.5 cross_entropy vs torch\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#@title Solution { display-mode: \"form\" }\n",
    "# solution: redefines softmax and cross_entropy; the checks below re-verify the reference.\n",
    "def softmax(logits, axis=-1):\n",
    "    logits = np.asarray(logits, dtype=float)\n",
    "    shifted = logits - logits.max(axis=axis, keepdims=True)\n",
    "    e = np.exp(shifted)\n",
    "    return e / e.sum(axis=axis, keepdims=True)\n",
    "\n",
    "def cross_entropy(logits, targets):\n",
    "    logits = np.asarray(logits, dtype=float)\n",
    "    targets = np.asarray(targets)\n",
    "    probs = softmax(logits, axis=-1)\n",
    "    N = logits.shape[0]\n",
    "    true_class_probs = probs[np.arange(N), targets]\n",
    "    return float(-np.log(true_class_probs).mean())\n",
    "\n",
    "check(\"0.5 softmax distribution\", _softmax_is_a_distribution, required=True)\n",
    "check(\"0.5 softmax no overflow\",  _softmax_survives_overflow, required=True)\n",
    "check(\"0.5 cross_entropy vs torch\", _cross_entropy_matches_torch, required=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50f1d7a2",
   "metadata": {},
   "source": [
    "> **Interpretation.** Your from-scratch cross-entropy reproduces `F.cross_entropy` to floating-point tolerance, which means you now understand the loss every classifier in this curriculum minimizes. The stable softmax is the reason it does not blow up on real logits.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76f08698",
   "metadata": {},
   "source": [
    "### KL divergence is cross-entropy minus entropy\n",
    "\n",
    "The Kullback-Leibler divergence $D_{KL}(p\\|q) = \\sum_x p(x)\\log\\frac{p(x)}{q(x)} = H(p,q) - H(p)$ measures how far $q$ is from $p$. It is non-negative and zero exactly when $p = q$, but it is not symmetric and not a distance. We verify both the identity and the zero-at-equality property.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "a8877a4f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.660336Z",
     "iopub.status.busy": "2026-06-10T18:44:02.660258Z",
     "iopub.status.idle": "2026-06-10T18:44:02.663664Z",
     "shell.execute_reply": "2026-06-10T18:44:02.663308Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "KL(p||q)        = 0.0405\n",
      "H(p,q) - H(p)   = 0.0405\n",
      "KL(p||p)        = 0.0000   (zero, as it must be)\n",
      "[ ok ] KL = H(p,q) - H(p), and KL(p||p) = 0\n"
     ]
    }
   ],
   "source": [
    "def kl(p, q):\n",
    "    p, q = np.asarray(p, float), np.asarray(q, float)\n",
    "    mask = p > 0\n",
    "    return float((p[mask] * np.log(p[mask] / q[mask])).sum())\n",
    "\n",
    "def cross_entropy_pq(p, q):\n",
    "    p, q = np.asarray(p, float), np.asarray(q, float)\n",
    "    mask = p > 0\n",
    "    return float(-(p[mask] * np.log(q[mask])).sum())\n",
    "\n",
    "p = np.array([0.3, 0.5, 0.2])\n",
    "q = np.array([0.2, 0.5, 0.3])\n",
    "print(f\"KL(p||q)        = {kl(p, q):.4f}\")\n",
    "print(f\"H(p,q) - H(p)   = {cross_entropy_pq(p, q) - entropy(p):.4f}\")\n",
    "print(f\"KL(p||p)        = {kl(p, p):.4f}   (zero, as it must be)\")\n",
    "assert np.isclose(kl(p, q), cross_entropy_pq(p, q) - entropy(p)), \"KL = cross-entropy minus entropy\"\n",
    "assert abs(kl(p, p)) < 1e-12, \"KL of a distribution from itself is exactly 0\"\n",
    "print(\"[ ok ] KL = H(p,q) - H(p), and KL(p||p) = 0\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdc4e0eb",
   "metadata": {},
   "source": [
    "> **Caveat:** $D_{KL}(p\\|q) \\neq D_{KL}(q\\|p)$ in general, and KL does not satisfy the triangle inequality, so it is not a metric. It is the \"distance\" you will see in VAEs, in RLHF objectives, and in mech-interp logit-difference work, always with the asymmetry mattering.\n",
    "\n",
    "> **Key takeaways.** Entropy is average surprise; uniform maximizes it, certainty zeroes it. Cross-entropy equals NLL for one-hot targets and is the universal classification loss. Always subtract the max in softmax. KL is cross-entropy minus entropy: non-negative, zero at equality, asymmetric.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c18f7e95",
   "metadata": {},
   "source": [
    "## Part 6 — Putting it together: gradient descent from scratch\n",
    "\n",
    "> **Objectives.** Assemble the pieces into the smallest complete ML algorithm: gradient descent on linear regression that recovers known parameters from noisy synthetic data, with a loss curve and an experiment log.\n",
    "\n",
    "This is the consolidation cell. Every idea so far converges here: a parameterized function, a loss, its gradient (which we grad-checked the form of in Part 3), and the update $\\mathbf{w} \\leftarrow \\mathbf{w} - \\eta\\nabla L$. We synthesize data from a known law `y = 3x + 5 + noise` so we can assert that the fit recovers `(3, 5)`: ground-truth recovery is the strongest check available for a learning algorithm.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "473791a8",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.664459Z",
     "iopub.status.busy": "2026-06-10T18:44:02.664393Z",
     "iopub.status.idle": "2026-06-10T18:44:02.667407Z",
     "shell.execute_reply": "2026-06-10T18:44:02.667017Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data shapes: (200, 1) (200,)   true params: w=3.0, b=5.0\n"
     ]
    }
   ],
   "source": [
    "def make_regression_data(n=200, true_w=3.0, true_b=5.0, noise=0.1, seed=SEED):\n",
    "    gen = np.random.default_rng(seed)\n",
    "    X = gen.standard_normal((n, 1))\n",
    "    y = true_w * X[:, 0] + true_b + noise * gen.standard_normal(n)\n",
    "    return X, y\n",
    "\n",
    "def mse_loss(w, b, X, y):\n",
    "    residual = X @ w + b - y           # (N,)\n",
    "    return float(0.5 * (residual ** 2).mean())\n",
    "\n",
    "def mse_gradient(w, b, X, y):\n",
    "    N = X.shape[0]\n",
    "    residual = X @ w + b - y           # (N,)\n",
    "    dw = X.T @ residual / N            # (d,) the X^T(Xw - y)/N from Part 3\n",
    "    db = float(residual.mean())        # bias adds to every prediction equally\n",
    "    return dw, db\n",
    "\n",
    "Xreg, yreg = make_regression_data()\n",
    "print(\"data shapes:\", Xreg.shape, yreg.shape, \"  true params: w=3.0, b=5.0\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "819ae436",
   "metadata": {},
   "source": [
    "The training loop carries the four-comment skeleton you will see in every deep-learning chapter: forward, backward, update, track. We watch the convergence logic honestly: the loop runs a fixed `STEPS` count with no early break, so it cannot exit on a spurious first-iteration condition. We log the loss every so often as the expected-value reference.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "aa4a5d72",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.668308Z",
     "iopub.status.busy": "2026-06-10T18:44:02.668236Z",
     "iopub.status.idle": "2026-06-10T18:44:02.676374Z",
     "shell.execute_reply": "2026-06-10T18:44:02.676075Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "recovered w=2.9931, b=4.9913   (true 3.0, 5.0)\n",
      "loss: start 16.8288 -> end 0.005234   over 500 steps\n",
      "[ ok ] gradient descent recovered the ground-truth parameters\n"
     ]
    }
   ],
   "source": [
    "def gradient_descent(X, y, lr=0.1, n_steps=STEPS):\n",
    "    w = np.zeros(X.shape[1]); b = 0.0\n",
    "    losses = []\n",
    "    for step in range(n_steps):\n",
    "        # track\n",
    "        losses.append(mse_loss(w, b, X, y))\n",
    "        # backward (no separate forward: the loss/grad recompute the prediction)\n",
    "        dw, db = mse_gradient(w, b, X, y)\n",
    "        # update\n",
    "        w = w - lr * dw\n",
    "        b = b - lr * db\n",
    "    return w, b, losses\n",
    "\n",
    "w_hat, b_hat, losses = gradient_descent(Xreg, yreg, lr=0.1, n_steps=STEPS)\n",
    "print(f\"recovered w={w_hat[0]:.4f}, b={b_hat:.4f}   (true 3.0, 5.0)\")\n",
    "print(f\"loss: start {losses[0]:.4f} -> end {losses[-1]:.6f}   over {STEPS} steps\")\n",
    "assert abs(w_hat[0] - 3.0) < 0.1 and abs(b_hat - 5.0) < 0.1, \\\n",
    "    \"GD should recover the true (3, 5) within 0.1; if not, check the gradient or lr\"\n",
    "print(\"[ ok ] gradient descent recovered the ground-truth parameters\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "3991f7d9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.677158Z",
     "iopub.status.busy": "2026-06-10T18:44:02.677037Z",
     "iopub.status.idle": "2026-06-10T18:44:02.961623Z",
     "shell.execute_reply": "2026-06-10T18:44:02.961325Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAGGCAYAAACNCg6xAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAU5pJREFUeJzt3XlcVOX+B/DPDMuArCKbKIuKhbiAIRKmqYkSeU2tbmZWuJSVeC0pvXL7pXm1LLuZWlxNTU3LUtRc08TdyuuOmriLGyCgKJvKNs/vD5yT44w4gzNzBvi8Xy9eMOc8nPOd7wzw5Xme8xyFEEKAiIiIiB5IKXcARERERLUFCyciIiIiA7FwIiIiIjIQCyciIiIiA7FwIiIiIjIQCyciIiIiA7FwIiIiIjIQCyciIiIiA7FwIiIiIjIQCyeyiMGDByMoKEhrm0KhwEcffSRLPMaqTbESVWfq1KkICQmBWq2WO5R6KygoCIMHD5Y7DJNIT0+Hra0t/vzzT7lDsRgWTnVcRkYGRo4ciUceeQQNGjRAgwYNEBoaioSEBBw5ckTu8MxuyZIlmD59utxhWKVffvmFxWA9U1hYiM8++wz//Oc/oVT+9es/KCgICoVC5+Ott97SOcaNGzcwfPhweHl5wcnJCd27d8fBgwf1nm/NmjV47LHH4ODggICAAEyYMAEVFRUmeS58/z6YMa/Vvfbu3YsRI0YgIiICdnZ2UCgUetuFhoaid+/eGD9+vClDt2q2cgdA5rNu3ToMGDAAtra2GDRoEMLCwqBUKnHixAmsXLkSs2bNQkZGBgIDA2WJ79atW7C1Ne9bcMmSJfjzzz/x7rvvmvU8tdEvv/yC5ORk/vGpR+bPn4+KigoMHDhQZ194eDjee+89rW2PPPKI1mO1Wo3evXvj8OHDGDNmDDw9PfHf//4X3bp1w4EDB9CyZUup7YYNG9CvXz9069YNX331FY4ePYrJkycjNzcXs2bNeujnwvdv9Yx5rfT55ZdfMG/ePLRr1w7NmzfHqVOn7tv2rbfewjPPPIOzZ8+iRYsWpn4q1kdQnXTmzBnh5OQkWrVqJbKysnT2l5eXixkzZoiLFy9We5zi4mKTxBMfHy8CAwNNcixj9O7d2yTnBSAmTJjw0MexJgkJCYK/AoSorKwUt27dkjsMi2jXrp145ZVXdLYHBgaK3r17P/D7ly5dKgCIlJQUaVtubq5wd3cXAwcO1GobGhoqwsLCRHl5ubTtgw8+EAqFQhw/fvwhnkWV2vz+DQwMFPHx8dW2edjfvca8VvpcuXJF3Lx5Uwjx4FyXlZWJhg0big8//PChYq4tOFRXR02dOhUlJSVYsGABGjdurLPf1tYWo0aNgr+/v7Rt8ODBcHZ2xtmzZ/HMM8/AxcUFgwYNAgDs2rULf//73xEQEACVSgV/f3+MHj0at27d0jn2qlWr0KZNGzg4OKBNmzb4+eef9caob95QZmYmhg4dCh8fH6hUKrRu3Rrz58/XarN9+3YoFAosW7YMH3/8MZo2bQoHBwf06NEDZ86ckdp169YN69evx4ULF6Shh3vnWd2rtLQUo0ePhpeXF1xcXPDss8/i8uXLetsaEisAfPXVV2jdujUaNGiAhg0bokOHDliyZInOsYYNGwY/Pz+oVCo0a9YMb7/9NsrKyqQ2N27cwLvvvgt/f3+oVCoEBwfjs88+05qrcv78eSgUCvznP//BnDlz0KJFC6hUKkRGRmLfvn1Su8GDByM5OVl6HTQfD7JhwwZ07doVLi4ucHV1RWRkpM5zSUlJQUREBBwdHeHp6YlXXnkFmZmZWm0077XMzEz069cPzs7O8PLywvvvv4/KykoAQHl5OTw8PDBkyBCdOAoLC+Hg4ID3339f2lZaWooJEyYgODhYeo+OHTsWpaWlWt+rUCgwcuRI/PDDD2jdujVUKhU2btwIADhy5Ai6du0KR0dHNG3aFJMnT8aCBQugUChw/vx5nVx06dIFTk5OcHFxQe/evXHs2DGjn6eGWq3GjBkz0LZtWzg4OMDLywtPP/009u/fr9Xu+++/l/Lr4eGBl156CZcuXbrfSybJyMjAkSNHEBMTc982ZWVlKCkpue/+5cuXw8fHB88995y0zcvLCy+++CJWr14t5To9PR3p6ekYPny4Vq/yiBEjIITA8uXLq421vLwcEydORMuWLeHg4IBGjRqhc+fOSE1NBfDg969arcb06dPRunVrODg4wMfHB2+++SauX7+udZ6goCD87W9/w6ZNmxAeHg4HBweEhoZi5cqVRsXzsBYuXAiFQoEdO3ZgxIgR8Pb2RtOmTR/qmIa+Vvfj4+MDR0dHg85lZ2eHbt26YfXq1Q8Vc23Bobo6at26dQgODkZUVJRR31dRUYHY2Fh07twZ//nPf9CgQQMAVX8Mb968ibfffhuNGjXC3r178dVXX+Hy5ctISUmRvn/Tpk14/vnnERoaiilTpuDatWsYMmSIQb8EcnJy8Pjjj0t/2Ly8vLBhwwYMGzYMhYWFOsNtn376KZRKJd5//30UFBRg6tSpGDRoEPbs2QMA+OCDD1BQUIDLly/jyy+/BAA4OztXG8Prr7+O77//Hi+//DI6deqErVu3onfv3jWOde7cuRg1ahReeOEFvPPOO7h9+zaOHDmCPXv24OWXXwYAZGVloWPHjtJ8hJCQEGRmZmL58uW4efMm7O3tcfPmTXTt2hWZmZl48803ERAQgD/++ANJSUnIzs7Wmce1ZMkSFBUV4c0334RCocDUqVPx3HPP4dy5c7Czs8Obb76JrKwspKamYvHixQ98bYCqX+5Dhw5F69atkZSUBHd3dxw6dAgbN26UnsvChQsxZMgQREZGYsqUKcjJycGMGTPw+++/49ChQ3B3d5eOV1lZidjYWERFReE///kPNm/ejC+++AItWrTA22+/DTs7O/Tv3x8rV67EN998A3t7e+l7V61ahdLSUrz00ksAqv5QPvvss/jtt98wfPhwtGrVCkePHsWXX36JU6dOYdWqVVrPZevWrVi2bBlGjhwJT09PBAUFITMzE927d4dCoUBSUhKcnJwwb948qFQqnVwsXrwY8fHxiI2NxWeffYabN29i1qxZ6Ny5Mw4dOqRVoD/oeWoMGzYMCxcuRFxcHF5//XVUVFRg165d+N///ocOHToAAD7++GN8+OGHePHFF/H6668jLy8PX331FZ588kmd/N7rjz/+AAA89thjevdv3boVDRo0QGVlJQIDAzF69Gi88847Wm0OHTqExx57TGt+FAB07NgRc+bMwalTp9C2bVscOnQIAKS4Nfz8/NC0aVNp//189NFHmDJlCl5//XV07NgRhYWF2L9/Pw4ePIiePXs+8P375ptvSu/FUaNGISMjA19//TUOHTqE33//HXZ2dlLb06dPY8CAAXjrrbcQHx+PBQsW4O9//zs2btyInj17GhSPqYwYMQJeXl4YP368VMCWl5ejoKDAoO/38PCQXhtDXytTiYiIwOrVq1FYWAhXV1eTHdcqyd3lRaZXUFAgAIh+/frp7Lt+/brIy8uTPjRdsUJUDacBEOPGjdP5vrvbaUyZMkUoFApx4cIFaVt4eLho3LixuHHjhrRt06ZNAoDOkBnuGf4aNmyYaNy4sbh69apWu5deekm4ublJMWzbtk0AEK1atRKlpaVSuxkzZggA4ujRo9I2Y4bq0tLSBAAxYsQIre0vv/xyjWPt27evaN26dbXnfe2114RSqRT79u3T2adWq4UQQkyaNEk4OTmJU6dOae0fN26csLGxkYZcMzIyBADRqFEjkZ+fL7VbvXq1ACDWrl0rbTNmqOPGjRvCxcVFREVF6QxraWIsKysT3t7eok2bNlpt1q1bJwCI8ePHS9s077V///vfWsdq3769iIiIkB7/+uuvOnELIcQzzzwjmjdvLj1evHixUCqVYteuXVrtZs+eLQCI33//XdoGQCiVSnHs2DGttv/4xz+EQqEQhw4dkrZdu3ZNeHh4CAAiIyNDCCFEUVGRcHd3F2+88YbW91+5ckW4ublpbTf0eW7dulUAEKNGjRL30uT3/PnzwsbGRnz88cda+48ePSpsbW11tt/r//7v/wQAUVRUpLOvT58+4rPPPhOrVq0S3377rejSpYsAIMaOHavVzsnJSQwdOlTn+9evXy8AiI0bNwohhPj8888FAL1TASIjI8Xjjz9ebaxhYWEPHDq83/t3165dAoD44YcftLZv3LhRZ3tgYKAAIFasWCFtKygoEI0bNxbt27c3Kh5j3DtUt2DBAgFAdO7cWVRUVGi11fy+M+RD8x4VwvDXyhCG/K5YsmSJACD27Nlj8HFrKw7V1UGFhYUA9PeudOvWDV5eXtKHprv7bnf/F6xxd5dtSUkJrl69ik6dOkEIIf33mJ2djbS0NMTHx8PNzU1q37NnT4SGhlYbsxACK1asQJ8+fSCEwNWrV6WP2NhYFBQU6FwNMmTIEK1eiC5dugAAzp07V+257ueXX34BAIwaNUpr+709XcbE6u7ujsuXL2sNk91NrVZj1apV6NOnj85/5wCk4YeUlBR06dIFDRs21DpfTEwMKisrsXPnTq3vGzBgABo2bCg9ftjcpKamoqioCOPGjYODg4PeGPfv34/c3FyMGDFCq03v3r0REhKC9evX6xz33qu2unTpohXjU089BU9PTyxdulTadv36daSmpmLAgAHStpSUFLRq1QohISFa+XnqqacAANu2bdM6T9euXXXekxs3bkR0dDTCw8OlbR4eHtJw9d25uHHjBgYOHKh1LhsbG0RFRemcy5DnuWLFCigUCkyYMEHnezX5XblyJdRqNV588UWt8/r6+qJly5Z6z3u3a9euwdbWVu/vhTVr1mDs2LHo27cvhg4dih07diA2NhbTpk3TGqq+deuW3h44zeutGbrXfL5fW31D/Hdzd3fHsWPHcPr06Wrb6ZOSkgI3Nzf07NlTK08RERFwdnbWyZOfnx/69+8vPXZ1dcVrr72GQ4cO4cqVKw8djzHeeOMN2NjYaG0LCwtDamqqQR++vr7S9xn6WpmK5vfN1atXTXpca8ShujrIxcUFAFBcXKyz75tvvkFRURFycnLwyiuv6Oy3tbXVO6x28eJFjB8/HmvWrNGZJ6DpRr5w4QIA6L1a49FHH632Mti8vDzcuHEDc+bMwZw5c/S2yc3N1XocEBCg9Vjzg3tvfIa6cOEClEqlzlUhjz76aI1j/ec//4nNmzejY8eOCA4ORq9evfDyyy/jiSeekI5VWFiINm3aVBvb6dOnceTIEXh5eVV7Pg1T5+bs2bMAUG2cmtf/3nwBQEhICH777TetbZp5PPfGeXeMtra2eP7557FkyRKUlpZCpVJh5cqVKC8v1yqcTp8+jePHjxucn2bNmumNPzo6Wmd7cHCw1mPNH09NUXave4cpDHmeZ8+ehZ+fHzw8PPQeU3NeIcR9r4a6e/jpYSkUCowePRq//vortm/fLv2ucHR01Ds35vbt29L+uz/fr+2D5s78+9//Rt++ffHII4+gTZs2ePrpp/Hqq6+iXbt2D4z99OnTKCgogLe3t979974XgoODdeb3aa4mPH/+PHx9fR8qHmPoe182bNiw2nlp92Poa2UqQggAMGiuZG3HwqkOcnNzQ+PGjfUuSKaZ83TvRFcNlUqlMyZeWVmJnj17Ij8/H//85z8REhICJycnZGZmYvDgwSZZSE9zjFdeeQXx8fF629z7S+re/8w0ND/A5mJMrK1atcLJkyexbt06bNy4EStWrMB///tfjB8/HhMnTjTqnD179sTYsWP17r/3snG5cmOM+8V4r5deegnffPONdHn7smXLEBISgrCwMKmNWq1G27ZtMW3aNL3HuPsiCODh/mhoXv/Fixdr/Yevce8SG4Y+T0POq1AosGHDBr3HfND8vUaNGqGiogJFRUXSP1fV0eQsPz9f2ta4cWNkZ2frtNVs8/Pzk9pptt+b++zsbHTs2LHacz/55JM4e/YsVq9ejU2bNmHevHn48ssvMXv2bLz++uvVfq9arYa3tzd++OEHvfvvV1ybKx5j6HtflpWVab0G1fHy8pLeG4a+Vqai+UfA09PTpMe1Riyc6qjevXtj3rx52Lt37wN/ST3I0aNHcerUKXz33Xd47bXXpO33XlGiWQ9KX3f2yZMnqz2H5iq2ysrKGv13dT/G/PcTGBgItVqNs2fPavWa3Bu7sbE6OTlhwIABGDBgAMrKyvDcc8/h448/RlJSEry8vODq6vrAVXdbtGiB4uJi2XKj6YX7888/dXpgNDSv/8mTJ3V6Y06ePFnj9cKefPJJNG7cGEuXLkXnzp2xdetWfPDBBzrxHT58GD169Kjxf7yBgYFaV2Vq3LtNkwtvb2+TvR4tWrTAr7/+ivz8/Pv2OrVo0QJCCDRr1kynUDZESEgIgKqr6wzpKdEMJd5daISHh2PXrl1Qq9Va/2Dt2bMHDRo0kOLSDHfu379f6/dPVlYWLl++jOHDhz/w/JorKocMGYLi4mI8+eST+Oijj6RC5X6vc4sWLbB582Y88cQTBhXIZ86cgRBC63iaNYvunuT/oHjM5Y8//kD37t0NapuRkSHFbOhrZSoZGRlQKpUmP6414hynOmrs2LFo0KABhg4dipycHJ39xvQ8aP6Duft7hBCYMWOGVrvGjRsjPDwc3333ndZVIKmpqUhPT3/gOZ5//nmsWLFCbxGRl5dncLx3c3JyMviKlLi4OADAzJkztbbfe8WaMbFeu3ZNa5+9vT1CQ0MhhEB5eTmUSiX69euHtWvX6lx2DvyV8xdffBG7d+/Gr7/+qtPmxo0bNVqN2cnJSfr+B+nVqxdcXFwwZcoUqav/3hg7dOgAb29vzJ49W2uIYMOGDTh+/LjeqxMNoVQq8cILL2Dt2rVYvHgxKioqtIbpgKr8ZGZmYu7cuTrff+vWrWovsdeIjY3F7t27kZaWJm3Lz8/X6bmIjY2Fq6srPvnkE5SXl+scpybv1eeffx5CCL29kJr8Pvfcc7CxscHEiRN1fn6FEDrvtXtphiHvfZ/l5+frLI1QXl6OTz/9FPb29lp/tF944QXk5ORoXa5/9epVpKSkoE+fPtKcmtatWyMkJARz5szROvasWbOgUCjwwgsvVBvrvc/F2dkZwcHBWu+r+71/X3zxRVRWVmLSpEk6x62oqNBpn5WVpbVkSmFhIRYtWoTw8HCpR9GQeMylpnOcDH2tgKqhYs1wfE0dOHAArVu31prfWlexx6mOatmyJZYsWYKBAwfi0UcflVYOF0IgIyMDS5YsgVKpNGiZgJCQELRo0QLvv/8+MjMz4erqihUrVuidLzNlyhT07t0bnTt3xtChQ5Gfny+tY6RvztXdPv30U2zbtg1RUVF44403EBoaivz8fBw8eBCbN282uLv6bhEREVi6dCkSExMRGRkJZ2dn9OnTR2/b8PBwDBw4EP/9739RUFCATp06YcuWLXp7IQyNtVevXvD19cUTTzwBHx8fHD9+HF9//TV69+4tDZd88skn2LRpE7p27SpdSp+dnY2UlBT89ttvcHd3x5gxY7BmzRr87W9/w+DBgxEREYGSkhIcPXoUy5cvx/nz543uIo+IiABQNRk+NjYWNjY20uX993J1dcWXX36J119/HZGRkXj55ZfRsGFDHD58GDdv3sR3330HOzs7fPbZZxgyZAi6du2KgQMHSssRBAUFYfTo0UbFd7cBAwbgq6++woQJE9C2bVu0atVKa/+rr76KZcuW4a233sK2bdvwxBNPoLKyEidOnMCyZcvw66+/6p18f7exY8fi+++/R8+ePfGPf/xDWo4gICAA+fn5Uo+Eq6srZs2ahVdffRWPPfYYXnrpJXh5eeHixYtYv349nnjiCXz99ddGPb/u3bvj1VdfxcyZM3H69Gk8/fTTUKvV2LVrF7p3746RI0eiRYsWmDx5MpKSknD+/Hn069cPLi4uyMjIwM8//4zhw4drrWt1r+bNm6NNmzbYvHkzhg4dKm1fs2YNJk+ejBdeeAHNmjVDfn6+tOL+J598ovPH+PHHH8eQIUOQnp4urUZdWVmpU/R9/vnnePbZZ9GrVy+89NJL+PPPP/H111/j9ddf13n97hUaGopu3bohIiICHh4e2L9/P5YvX46RI0dKbe73/u3atSvefPNNTJkyBWlpaejVqxfs7Oxw+vRppKSkYMaMGVqF2yOPPIJhw4Zh37598PHxwfz585GTk4MFCxYYFc/58+fRrFkzxMfHY+HChdU+P2PUdI6TMa9Vjx49AGhP4bhw4YK01IOm2J48eTKAqt7ZV199VWpbXl4urUFVL1jyEj6yvDNnzoi3335bBAcHCwcHB+Ho6ChCQkLEW2+9JdLS0rTaxsfHCycnJ73HSU9PFzExMcLZ2Vl4enqKN954Qxw+fFgAEAsWLNBqu2LFCtGqVSuhUqlEaGioWLlypd6Vw6FnNe6cnByRkJAg/P39hZ2dnfD19RU9evQQc+bMkdpoLs+9e0VcIf66FP/ueIqLi8XLL78s3N3d9S6JcK9bt26JUaNGiUaNGgknJyfRp08fcenSpRrH+s0334gnn3xSNGrUSKhUKtGiRQsxZswYUVBQoHWsCxcuiNdee014eXkJlUolmjdvLhISErSWWygqKhJJSUkiODhY2NvbC09PT9GpUyfxn//8R5SVlWnl4PPPP9d5bvc+h4qKCvGPf/xDeHl5CYVCYdDSBGvWrBGdOnUSjo6OwtXVVXTs2FH8+OOPWm2WLl0q2rdvL1QqlfDw8BCDBg0Sly9f1mpzv/fahAkT9MahVquFv7+/ACAmT56sN7aysjLx2WefidatWwuVSiUaNmwoIiIixMSJE7XyDUAkJCToPcahQ4dEly5dhEqlEk2bNhVTpkwRM2fOFADElStXtNpu27ZNxMbGCjc3N+Hg4CBatGghBg8eLPbv31+j51lRUSE+//xzERISIuzt7YWXl5eIi4sTBw4c0Gq3YsUK0blzZ+Hk5CScnJxESEiISEhIECdPntT7nO42bdo04ezsrLW8yP79+0WfPn1EkyZNhL29vXB2dhadO3cWy5Yt03uM/Px8MWzYMNGoUSPRoEED0bVrV71LaQghxM8//yzCw8OlfP7f//2f9F6tzuTJk0XHjh2Fu7u79Dvr448/1vreB71/58yZIyIiIoSjo6NwcXERbdu2FWPHjtW6k4JmxfRff/1VtGvXTqhUKhESEqLzu8WQeI4ePXrf5Vzudb/lCO6Xx5oy9LUKDAzU+d1Y3TIIXbt21Wq7YcMGAUCcPn3apPFbK4UQVjRblIjIyrz77rv45ptvUFxcbLKJ3nIpKChA8+bNMXXqVAwbNkzucGQXFBSENm3aYN26dQ99rP/+978YO3Yszp49Cx8fHxNEV3v069cPCoXivneJqGs4x4mI6I5717a5du0aFi9ejM6dO9f6ogmouuJ27Nix+Pzzz01yNSz9Zdu2bRg1alS9K5qOHz+OdevW6Z1TVlexx4mI6I7w8HB069YNrVq1Qk5ODr799ltkZWVhy5YtePLJJ+UOj0zMlD1OVH9wcjgR0R3PPPMMli9fjjlz5kChUOCxxx7Dt99+y6KJiCTscSIiIiIyEOc4ERERERmIhRMRERGRgTjHqRpqtRpZWVlwcXGpFzcuJCIiqo+EECgqKoKfn5/O/VrvxcKpGllZWTo3qCQiIqK66dKlSw+8owYLp2pobolx6dIluLq6mvz4arUaeXl58PLyemCFSw+P+bYs5tvymHPLYr4tz1w5LywshL+/v/R3vzosnKpx972pzFU43b59G66urvyhswDm27KYb8tjzi2L+bY8c+fckGk5fKWJiIiIDMTCiYiIiMhALJyIiIiIDMTCiYiIiMhALJyIiIiIDMTCSY/k5GSEhoYiMjJS7lCIiIjIirBw0iMhIQHp6enYt2+f3KEQERGRFWHhRERERGQgFk5EREREBmLhRERERGQgFk4ymf0T8MRABRb87CR3KERERGQgFk4yKa8AsnIVuHSFtwskIiKqLVg4ycTft+pzVq6NvIEQERGRwVg4ycS/cdXnrDwWTkRERLUFCyeZBNwpnHLzlSgrlzcWIiIiMgwLJ5l4NgQcHQSEUCArV+5oiIiIyBAsnGSiUABNfaq+vpgtbyxERERkGBZOMtJMEL/MwomIiKhWYOEko6Z35jlduqKQNxAiIiIyCAsnGQU0FgCAS1dkDoSIiIgMwsJJRpqhukscqiMiIqoVWDjJSBqqY+FERERUK7BwkpGmx+l6oQLFN+WNhYiIiB6sXhRO/fv3R8OGDfHCCy/IHYoWFyfA1VkNgL1OREREtUG9KJzeeecdLFq0SO4w9GriXQmAhRMREVFtUC8Kp27dusHFxUXuMPTy86oAwCvriIiIagOrL5x27tyJPn36wM/PDwqFAqtWrdJpk5ycjKCgIDg4OCAqKgp79+61fKA15Henx+lilsyBEBER0QNZfeFUUlKCsLAwJCcn692/dOlSJCYmYsKECTh48CDCwsIQGxuL3NzacQM4P687Q3XscSIiIrJ6tnIH8CBxcXGIi4u77/5p06bhjTfewJAhQwAAs2fPxvr16zF//nyMGzfOqHOVlpaitLRUelxYWAgAUKvVUKvVNYi+emq1Go01Q3XZAmq1MPk56C9qtRpCCLO8lqSL+bY85tyymG/LM1fOjTme1RdO1SkrK8OBAweQlJQkbVMqlYiJicHu3buNPt6UKVMwceJEne15eXm4ffv2Q8Wqj1qthovjTQCNcDFbICcnFwrefcVs1Go1CgoKIISAUmn1na21HvNtecy5ZTHflmeunBcVFRnctlYXTlevXkVlZSV8fHy0tvv4+ODEiRPS45iYGBw+fBglJSVo2rQpUlJSEB0drXO8pKQkJCYmSo8LCwvh7+8PLy8vuLq6mjx+tVqN8oo8KBQCt0uVsLH3hmdDk5+G7lCr1VAoFPDy8uIvOQtgvi2PObcs5tvyzJVzBwcHg9vW6sLJUJs3bzaonUqlgkql0tmuVCrN9kOhslfA1xPIzgMu5yjh3cgsp6E7FAqFWV9P0sZ8Wx5zblnMt+WZI+fGHKtWv9Kenp6wsbFBTk6O1vacnBz4+vrKFJXxmvKedURERLVCrS6c7O3tERERgS1btkjb1Go1tmzZoncozlDJyckIDQ1FZGSkKcJ8IN7sl4iIqHaw+qG64uJinDlzRnqckZGBtLQ0eHh4ICAgAImJiYiPj0eHDh3QsWNHTJ8+HSUlJdJVdjWRkJCAhIQEFBYWws3NzRRPo1oBfgKAgksSEBERWTmrL5z279+P7t27S481k7fj4+OxcOFCDBgwAHl5eRg/fjyuXLmC8PBwbNy4UWfCuDXjUB0REVHtYPWFU7du3SBE9esbjRw5EiNHjrRQRKanGaq7yMKJiIjIqtXqOU7mYvE5To2rPmflABWVFjklERER1QALJz0SEhKQnp6Offv2WeR8Po0AezugUl21LAERERFZJxZOVkCp5DwnIiKi2oCFk5WQ5jllyRsHERER3R8LJz0sPccJ+GueE5ckICIisl4snPSw9BwngItgEhER1QYsnKwEe5yIiIisHwsnKxGgKZzY40RERGS1WDhZCU2PU14+cOu2vLEQERGRfiycrISbC+DiVPU1h+uIiIisEwsnPeS4qk6h4ARxIiIia8fCSQ85rqoDOEGciIjI2rFwsiL+nCBORERk1Vg4WREWTkRERNaNhZMV4RwnIiIi68bCyYrcPcdJCHljISIiIl0snPSQ46o64K8ep6ISoKDIoqcmIiIiA7Bw0kOuq+ocHQAvj6qvOVxHRERkfVg4WRnNcN1FFk5ERERWh4WTleEEcSIiIuvFwsnKcBFMIiIi68XCycqwx4mIiMh6sXCyMgF+VZ/Z40RERGR9WDjpIddyBMBfQ3WXrwBqtcVPT0RERNVg4aSHXMsRAEBjL8BGCZSVAznXLH56IiIiqgYLJytjawP4+VR9zXlORERE1oWFkxXiBHEiIiLrxMLJCnGCOBERkXVi4WSFND1OF7PkjYOIiIi0sXCyQlwEk4iIyDqxcLJCUuHEOU5ERERWhYWTFQq4UzhduQqUlskbCxEREf2FhZMVauQOODoAQgBZuXJHQ0RERBosnPSQc+VwAFAoOEGciIjIGrFw0kPOlcM1OEGciIjI+rBwslJcBJOIiMj6sHCyUtIimCyciIiIrAYLJysl9ThxqI6IiMhqsHCyUpo5ThfZ40RERGQ1WDhZKU3hdKMQKCqRNxYiIiKqwsLJSjk3ABq6Vn3NeU5ERETWgYWTFZMmiHOeExERkVVg4WTFuCQBERGRdWHhZMU4QZyIiMi6sHCyYuxxIiIisi4snKyYdNsVFk5ERERWwdaYxsePH8dPP/2EXbt24cKFC7h58ya8vLzQvn17xMbG4vnnn4dKpTJXrPXO3ZPDhai6+S8RERHJx6Aep4MHDyImJgbt27fHb7/9hqioKLz77ruYNGkSXnnlFQgh8MEHH8DPzw+fffYZSktLzR23WSUnJyM0NBSRkZGyxuHnXVUs3S4F8q7LGgoRERHBwB6n559/HmPGjMHy5cvh7u5+33a7d+/GjBkz8MUXX+Bf//qXqWK0uISEBCQkJKCwsBBubm6yxWFvBzT2ArJyq4brvD1kC4WIiIhgYOF06tQp2NnZPbBddHQ0oqOjUV5e/tCBURV/378Kp4jWckdDRERUvxk0VHe/oun27dtGtSfjSRPEuQgmERGR7Iy+qk6tVmPSpElo0qQJnJ2dce7cOQDAhx9+iG+//dbkAdZ3AbyyjoiIyGoYXThNnjwZCxcuxNSpU2Fvby9tb9OmDebNm2fS4IiLYBIREVkTowunRYsWYc6cORg0aBBsbGyk7WFhYThx4oRJgyOu5URERGRNjC6cMjMzERwcrLNdrVZzUrgZaFYPz84FKirljYWIiKi+M7pwCg0Nxa5du3S2L1++HO3btzdJUPQX70aAyg6oVFddXUdERETyMWrlcAAYP3484uPjkZmZCbVajZUrV+LkyZNYtGgR1q1bZ44Y6zWlEmjqC5y9VDVcp5ksTkRERJZndI9T3759sXbtWmzevBlOTk4YP348jh8/jrVr16Jnz57miLHea8oJ4kRERFbB6B4nAOjSpQtSU1NNHQvdh2aeEyeIExERycvoHieyPF5ZR0REZB0M6nFq2LAhFAqFQQfMz89/qIBIF1cPJyIisg4GFU7Tp083cxhUHc2E8MvscSIiIpKVQYVTfHy8ueOgamgKp7zrwM1bQANHeeMhIiKqr2o0OVzj9u3bKCsr09rm6ur6UAGRLjcXwNUJKCwBLl8BHmkmd0RERET1k9GTw0tKSjBy5Eh4e3vDyckJDRs21Pog82jKeU5ERESyM7pwGjt2LLZu3YpZs2ZBpVJh3rx5mDhxIvz8/LBo0SJzxPhQ1q1bh0cffRQtW7as1Tch5s1+iYiI5Gf0UN3atWuxaNEidOvWDUOGDEGXLl0QHByMwMBA/PDDDxg0aJA54qyRiooKJCYmYtu2bXBzc0NERAT69++PRo0ayR2a0QK4JAEREZHsjO5xys/PR/PmzQFUzWfSLD/QuXNn7Ny507TRPaS9e/eidevWaNKkCZydnREXF4dNmzbJHVaNcBFMIiIi+RldODVv3hwZGRkAgJCQECxbtgxAVU+Uu7u7SYPbuXMn+vTpAz8/PygUCqxatUqnTXJyMoKCguDg4ICoqCjs3btX2peVlYUmTZpIj5s0aYLMzEyTxmgpXASTiIhIfkYXTkOGDMHhw4cBAOPGjUNycjIcHBwwevRojBkzxqTBlZSUICwsDMnJyXr3L126FImJiZgwYQIOHjyIsLAwxMbGIjc316RxWIO7F8EUQt5YiIiI6iuj5ziNHj1a+jomJgYnTpzAgQMHEBwcjHbt2pk0uLi4OMTFxd13/7Rp0/DGG29gyJAhAIDZs2dj/fr1mD9/PsaNGwc/Pz+tHqbMzEx07NjxvscrLS1FaWmp9LiwsBAAoFaroVarH/bp6FCr1RBCGHTsJt4AoETxTSD/hhoN3UweTp1nTL7p4THflsecWxbzbXnmyrkxx3uodZwAIDAwEIGBgQ97GKOVlZXhwIEDSEpKkrYplUrExMRg9+7dAICOHTvizz//RGZmJtzc3LBhwwZ8+OGH9z3mlClTMHHiRJ3teXl5uH37tsmfg1qtRkFBAYQQUCof3Pnn6e6FqzdscCQ9H61aVJg8nrrO2HzTw2G+LY85tyzm2/LMlfOioiKD2xpdOI0aNQrBwcEYNWqU1vavv/4aZ86csdjtWa5evYrKykr4+Phobffx8cGJEycAALa2tvjiiy/QvXt3qNVqjB07ttor6pKSkpCYmCg9LiwshL+/P7y8vMyysKdarYZCoYCXl5dBb4DAJgpcvQEUlXrA29vk4dR5xuabHg7zbXnMuWUx35Znrpw7ODgY3NbowmnFihVYs2aNzvZOnTrh008/tbr72j377LN49tlnDWqrUqmgUql0tiuVSrP9UCgUCoOP798YOHAMyMxRgj+jNWNMvunhMd+Wx5xbFvNteebIuTHHMvqs165dg5ub7gQbV1dXXL161djD1ZinpydsbGyQk5OjtT0nJwe+vr4Wi8OSeGUdERGRvIwunIKDg7Fx40ad7Rs2bJDWd7IEe3t7REREYMuWLdI2tVqNLVu2IDo6+qGOnZycjNDQUERGRj5smCbFtZyIiIjkZfRQXWJiIkaOHIm8vDw89dRTAIAtW7bgiy++MPkwXXFxMc6cOSM9zsjIQFpaGjw8PBAQEIDExETEx8ejQ4cO6NixI6ZPn46SkhLpKruaSkhIQEJCAgoLC/X2rsklwK/qM+9XR0REJA+jC6ehQ4eitLQUH3/8MSZNmgQACAoKwqxZs/Daa6+ZNLj9+/eje/fu0mPNxO34+HgsXLgQAwYMQF5eHsaPH48rV64gPDwcGzdu1JkwXldoepwuXwEqKwEbG3njISIiqm8UQtR8OcW8vDw4OjrC2dnZlDFZDU2PU0FBgdmuqsvNzYW3t7dBE9MqK4FHYoGKSmD3UsCPV9YZxdh808Nhvi2PObcs5tvyzJVzY/7eG33WW7du4ebNmwAALy8vXLt2DdOnT6+194DTx1rnONnYAH53OtM4z4mIiMjyjC6c+vbti0WLFgEAbty4gY4dO+KLL75A3759MWvWLJMHKIeEhASkp6dj3759coeigxPEiYiI5GN04XTw4EF06dIFALB8+XL4+vriwoULWLRoEWbOnGnyAElbwF33rCMiIiLLMrpwunnzJlxcXAAAmzZtwnPPPQelUonHH38cFy5cMHmApE2zltPFLHnjICIiqo9qtI7TqlWrcOnSJfz666/o1asXACA3N9csE6jlYK1znIC7FsFkjxMREZHFGV04jR8/Hu+//z6CgoIQFRUlLTa5adMmtG/f3uQByoFznIiIiEgfo9dxeuGFF9C5c2dkZ2cjLCxM2t6jRw/079/fpMGRLk2PU841oLQMUNnLGw8REVF9YnThBAC+vr4694Pr2LGjSQKi6jVyBxo4ADdvA5k5QHN/uSMiIiKqP7hiVy2jUNw1QZzDdURERBbFwkkPa54cDtw1QZyFExERkUWxcNLDmieHA5wgTkREJBcWTrUQlyQgIiKSh9GTw9esWaN3u0KhgIODA4KDg9GsWbOHDozuL4CLYBIREcnC6MKpX79+UCgUEEJobddsUygU6Ny5M1atWoWGDRuaLFD6C3uciIiI5GH0UF1qaioiIyORmpqKgoICFBQUIDU1FVFRUVi3bh127tyJa9eu4f333zdHvIS/CqeCIqCwWN5YiIiI6hOje5zeeecdzJkzB506dZK29ejRAw4ODhg+fDiOHTuG6dOnY+jQoSYNlP7i5Ah4uAH5BVW9Tq2D5Y6IiIiofjC6x+ns2bN670nn6uqKc+fOAQBatmyJq1evPnx0MrH25QgALklAREQkB6MLp4iICIwZMwZ5eXnStry8PIwdO1YqNE6fPg1//9q7pLW1L0cA3DVBnIUTERGRxRg9VPftt9+ib9++aNq0qVQcXbp0Cc2bN8fq1asBAMXFxfi///s/00ZKWtjjREREZHlGF06PPvoo0tPTsWnTJpw6dUra1rNnTyiVVR1Y/fr1M2mQpIuLYBIREVlejW7yq1Qq8fTTT+Ppp582dTxkIC5JQEREZHk1Wjl8x44d6NOnD4KDgxEcHIxnn30Wu3btMnVsVI27h+ruWVKLiIiIzMTowun7779HTEwMGjRogFGjRmHUqFFwcHBAjx49sGTJEnPESHr4eQNKJVBaBuTmyx0NERFR/WD0UN3HH3+MqVOnYvTo0dK2UaNGYdq0aZg0aRJefvllkwYoh+TkZCQnJ6OyslLuUO7L3g5o7AVk5lT1Ovk0kjsiIiKius/oHqdz586hT58+OtufffZZZGRkmCQoudWG5QgAThAnIiKyNKMLJ39/f2zZskVn++bNm2v12k21ESeIExERWZbRQ3XvvfceRo0ahbS0NOm2K7///jsWLlyIGTNmmDxAuj+u5URERGRZRhdOb7/9Nnx9ffHFF19g2bJlAIBWrVph6dKl6Nu3r8kDpPvz5+rhREREFlWjdZz69++P/v37mzoWMpJmjtNlFk5EREQWUaN1nMg6aO5Xl5UHlFfIGwsREVF9YFCPU8OGDaFQKAw6YH4+FxWyFC8PQGUHlJYD2blAgJ/cEREREdVtBhVO06dPN3MYVBNKJdC0MXD2YtU8JxZORERE5mVQ4RQfH2/uOKiG/O8UTryyjoiIyPwMmuNUUlJi1EGNbW9tkpOTERoaisjISLlDeSAugklERGQ5BhVOwcHB+PTTT5Gdff+/zkIIpKamIi4uDjNnzjRZgHKoLSuHA1wEk4iIyJIMGqrbvn07/vWvf+Gjjz5CWFgYOnToAD8/Pzg4OOD69etIT0/H7t27YWtri6SkJLz55pvmjpvuYI8TERGR5RhUOD366KNYsWIFLl68iJSUFOzatQt//PEHbt26BU9PT7Rv3x5z585FXFwcbGxszB0z3UUzIZw9TkREROZn1AKYAQEBeO+99/Dee++ZKx4ykmao7up1oOQW4OQobzxERER1GRfArOXcnAE3l6qvOVxHRERkXiyc6gDNCuIXs+SNg4iIqK5j4VQHBN6Z53SBhRMREZFZsXCqAzQTxC9yqI6IiMisWDjVAQHscSIiIrIIowunjRs34rfffpMeJycnIzw8HC+//DKuX79u0uDIMJo5TpwcTkREZF5GF05jxoxBYWEhAODo0aN477338MwzzyAjIwOJiYkmD5AeTDPH6VI2UFkpbyxERER1mVHrOAFARkYGQkNDAQArVqzA3/72N3zyySc4ePAgnnnmGZMHSA/W2AuwtQHKK4ArV4EmPnJHREREVDcZ3eNkb2+PmzdvAgA2b96MXr16AQA8PDykniiyLBsboOmdW69wgjgREZH5GF04de7cGYmJiZg0aRL27t2L3r17AwBOnTqFpk2bmjxAOSQnJyM0NBSRkZFyh2IwaUmCTHnjICIiqsuMLpy+/vpr2NraYvny5Zg1axaaNGkCANiwYQOefvppkwcoh4SEBKSnp2Pfvn1yh2Iwza1X2ONERERkPkbPcQoICMC6det0tn/55ZcmCYhqJrCqfuXq4URERGZkdI/TwYMHcfToUenx6tWr0a9fP/zrX/9CWVmZSYMjwwWwx4mIiMjsjC6c3nzzTZw6dQoAcO7cObz00kto0KABUlJSMHbsWJMHSIbhbVeIiIjMz+jC6dSpUwgPDwcApKSk4Mknn8SSJUuwcOFCrFixwtTxkYE0c5xuFAIFxfLGQkREVFcZXTgJIaBWqwFULUegWbvJ398fV69eNW10ZDDnBoBnw6qvOc+JiIjIPIwunDp06IDJkydj8eLF2LFjh7QcQUZGBnx8uPKinHhlHRERkXkZXThNnz4dBw8exMiRI/HBBx8gODgYALB8+XJ06tTJ5AGS4TTznNjjREREZB5GL0fQrl07ravqND7//HPY2NiYJCiqGenKOhZOREREZmF04aRx4MABHD9+HAAQGhqKxx57zGRBUc0E8Mo6IiIiszK6cMrNzcWAAQOwY8cOuLu7AwBu3LiB7t2746effoKXl5epYyQDSUN1nONERERkFkbPcfrHP/6B4uJiHDt2DPn5+cjPz8eff/6JwsJCjBo1yhwxkoE0Q3VZOUB5hbyxEBER1UVG9zht3LgRmzdvRqtWraRtoaGhSE5ORq9evUwaHBnHuxGgsgdKy6qKJ81tWIiIiMg0jO5xUqvVsLOz09luZ2cnre9E8lAquSQBERGRORldOD311FN45513kJX11wzkzMxMjB49Gj169DBpcGQ83nqFiIjIfIwunL7++msUFhYiKCgILVq0QIsWLdCsWTMUFhbiq6++MkeMZAQuSUBERGQ+Rs9x8vf3x8GDB7F582acOHECANCqVSvExMSYPDgyHnuciIiIzKdG6zgpFAr07NkTPXv2NHU8ZtG/f39s374dPXr0wPLly+UOx6wCuCQBERGR2RhUOM2cOdPgA1rjkgTvvPMOhg4diu+++07uUMzO/66hOiEAhULeeIiIiOoSgwqnL7/80qCDKRQKqyycunXrhu3bt8sdhkVo5jgV3wSuFwIebvLGQ0REVJcYNDk8IyPDoI9z584ZHcDOnTvRp08f+Pn5QaFQYNWqVTptkpOTERQUBAcHB0RFRWHv3r1Gn6e+cFABPp5VX1/IlDcWIiKiusboq+pMraSkBGFhYUhOTta7f+nSpUhMTMSECRNw8OBBhIWFITY2Frm5uVKb8PBwtGnTRufj7iUT6hPeeoWIiMg8anyTX1OJi4tDXFzcffdPmzYNb7zxBoYMGQIAmD17NtavX4/58+dj3LhxAIC0tDSTxFJaWorS0lLpcWFhIYCqRT/NsbinWq2GEMLkx/b3VWDvEQXOZ6rBNUn/Yq58k37Mt+Ux55bFfFueuXJuzPFkL5yqU1ZWhgMHDiApKUnaplQqERMTg927d5v8fFOmTMHEiRN1tufl5eH27dsmP59arUZBQQGEEFAqTdf518jNCYALTp27jdzcQpMdt7YzV75JP+bb8phzy2K+Lc9cOS8qKjK4rVUXTlevXkVlZSV8fHy0tvv4+EhrSBkiJiYGhw8fRklJCZo2bYqUlBRER0frtEtKSkJiYqL0uLCwEP7+/vDy8oKrq2vNn8h9qNVqKBQKeHl5mfQN0Cq46nPudUd4ezuY7Li1nbnyTfox35bHnFsW82155sq5g4PhfyutunAylc2bNxvUTqVSQaVS6WxXKpVm+6FQKBQmP35Q06rPl7IVUCq5HsHdzJFvuj/m2/KYc8tivi3PHDk35lgGt5w6dSpu3bolPf7999+15gMVFRVhxIgRBp/YEJ6enrCxsUFOTo7W9pycHPj6+pr0XHdLTk5GaGgoIiMjzXYOc9IsSXDlKnC7TN5YiIiI6hKDC6ekpCStMcC4uDhkZv51vfvNmzfxzTffmDQ4e3t7REREYMuWLdI2tVqNLVu26B1qM5WEhASkp6dj3759ZjuHOTVyB5wcqxbAvHxF7miIiIjqDoOH6oQQ1T6uqeLiYpw5c0Z6nJGRgbS0NHh4eCAgIACJiYmIj49Hhw4d0LFjR0yfPh0lJSXSVXakS6Go6nU6fq5qBfHgALkjIiIiqhtkn+O0f/9+dO/eXXqsmZwdHx+PhQsXYsCAAcjLy8P48eNx5coVhIeHY+PGjToTxklbgF9V4cSb/RIREZmO7IVTt27dHth7NXLkSIwcOdJCEVXNcUpOTkZlZaXFzmlq0s1+WTgRERGZjFGF07x58+Ds7AwAqKiowMKFC+HpWXV/D2PWQLB2CQkJSEhIQGFhIdzcaufN3jSrh7PHiYiIyHQMLpwCAgIwd+5c6bGvry8WL16s04asQ7M7SxKcvyxvHERERHWJwYXT+fPnzRgGmVpQk6rPF7OBikrA1kbeeIiIiOoCrthVR/l5Ayo7oLwCyMp5cHsiIiJ6MIMLp927d2PdunVa2xYtWoRmzZrB29sbw4cP11oQszar7QtgAoBSCfjfmeeUweE6IiIikzC4cPr3v/+NY8eOSY+PHj2KYcOGISYmBuPGjcPatWsxZcoUswRpabV9AUyNZneG6zIyq29HREREhjG4cEpLS0OPHj2kxz/99BOioqIwd+5cJCYmYubMmVi2bJlZgqSa0dyz7gILJyIiIpMwuHC6fv261qKTO3bsQFxcnPQ4MjISly5dMm109FCkHicO1REREZmEwYWTj48PMjIyAABlZWU4ePAgHn/8cWl/UVER7OzsTB8h1Zimx+k8e5yIiIhMwuDC6ZlnnsG4ceOwa9cuJCUloUGDBujSpYu0/8iRI2jRooVZgrS0ujA5HPhrLadLd5YkICIioodjcOE0adIk2NraomvXrpg7dy7mzp0Le3t7af/8+fPRq1cvswRpaXVlcrivJ6CyryqaLl+ROxoiIqLaz+AFMD09PbFz504UFBTA2dkZNjbaKyqmpKRIt2Mh66BUVt165dT5quE6zaKYREREVDNGL4Dp5uamUzQBgIeHh1YPFFkHzXAdJ4gTERE9PIN7nIYOHWpQu/nz59c4GDI9TS8T71lHRET08AwunBYuXIjAwEC0b98eQghzxkQmJPU48co6IiKih2Zw4fT222/jxx9/REZGBoYMGYJXXnkFHh4e5oyNTEBakoA9TkRERA/N4DlOycnJyM7OxtixY7F27Vr4+/vjxRdfxK+//lrneqDqynIEwF+LYF6+UnXDXyIiIqo5oyaHq1QqDBw4EKmpqUhPT0fr1q0xYsQIBAUFobi42FwxWlxdWY4AALwbAQ4qoFLNJQmIiIgeltFX1UnfqFRCoVBACIHKSq6uaK2Uyr8miPPKOiIioodjVOFUWlqKH3/8ET179sQjjzyCo0eP4uuvv8bFixe5hpMV45V1REREpmHw5PARI0bgp59+gr+/P4YOHYoff/wRnp6e5oyNTIRX1hEREZmGwYXT7NmzERAQgObNm2PHjh3YsWOH3nYrV640WXBkGpoepwssnIiIiB6KwYXTa6+9BoVCYc5YyEzY40RERGQaRi2AWV8kJycjOTm5zkx616zldPkKUFYO2NvJGw8REVFtVeOr6uqyurQcAQB4ewANHAC1GriULXc0REREtRcLp3pAobjryjoO1xEREdUYC6d6okVA1eczF+WNg4iIqDZj4VRPSIXTBXnjICIiqs1YONUTwXcKp7PscSIiIqoxFk71RHBg1eczF4E6dk9mIiIii2HhVE80a1o1SbygCLh6Xe5oiIiIaicWTvWEgwrwb1z1NYfriIiIaoaFUz0SzCvriIiIHgoLJz2Sk5MRGhqKyMhIuUMxKV5ZR0RE9HBYOOlR11YO1+CVdURERA+HhVM9wkUwiYiIHg4Lp3qk5Z0lCbJygZJb8sZCRERUG7FwqkfcXQHPhlVfn2OvExERkdFYONUzLfyrPp9m4URERGQ0Fk71jLSCOK+sIyIiMhoLp3qGV9YRERHVHAuneoZX1hEREdUcC6d6RjNUd/4yUF4hbyxERES1DQuneqaxF+DkCFRUAhmX5Y6GiIiodmHhVM8olcAjzaq+PnFO3liIiIhqGxZO9VCr5lWfWTgREREZh4VTPRTCwomIiKhGWDjpkZycjNDQUERGRsodilm0alH1mYUTERGRcVg46ZGQkID09HTs27dP7lDM4tE7PU6ZOUBBsbyxEBER1SYsnOohN2fAz7vq65PsdSIiIjIYC6d6ivOciIiIjMfCqZ7ilXVERETGY+FUT4XcmSB+nIUTERGRwVg41VPSUN1ZoLJS3liIiIhqCxZO9VQLf8DRAbh5GzjHW68QEREZhIVTPWVjA7RpWfX1kRPyxkJERFRbsHCqx9o9WvX5yEl54yAiIqotWDjVY20fqfp85JS8cRAREdUWLJzqsbCQqs/HTgMVnCBORET0QCyc6rGgJoCLE1BaBpzKkDsaIiIi68fCqR5TKoE2muE6znMiIiJ6IBZO9Vy7O4XTUc5zIiIieiAWTvWc5sq6Q8fljYOIiKg2YOFUz3VoU/X5+FmgsFjeWIiIiKwdC6d6ztcLCPQD1GrgwDG5oyEiIrJudb5wunTpErp164bQ0FC0a9cOKSkpcodkdTq2q/q894i8cRAREVm7Ol842draYvr06UhPT8emTZvw7rvvoqSkRO6wrEpUWNXnPYfljYOIiMja2codgLk1btwYjRs3BgD4+vrC09MT+fn5cHJykjky6xF1p8fpyEngdingoJI3HiIiImsle4/Tzp070adPH/j5+UGhUGDVqlU6bZKTkxEUFAQHBwdERUVh7969NTrXgQMHUFlZCX9//4eMum7xbwz4egLlFby6joiIqDqyF04lJSUICwtDcnKy3v1Lly5FYmIiJkyYgIMHDyIsLAyxsbHIzc2V2oSHh6NNmzY6H1lZWVKb/Px8vPbaa5gzZ47Zn1Nto1D8Nc+Jw3VERET3J/tQXVxcHOLi4u67f9q0aXjjjTcwZMgQAMDs2bOxfv16zJ8/H+PGjQMApKWlVXuO0tJS9OvXD+PGjUOnTp2qbVdaWio9LiwsBACo1Wqo1WpDn5LB1Go1hBBmObaxotoBa7YqsXOfwKhXhdzhmIU15bs+YL4tjzm3LObb8syVc2OOJ3vhVJ2ysjIcOHAASUlJ0jalUomYmBjs3r3boGMIITB48GA89dRTePXVV6ttO2XKFEycOFFne15eHm7fvm1c8AZQq9UoKCiAEAJKpbydf22DlQC8cTAdOHkmDw1d617xZE35rg+Yb8tjzi2L+bY8c+W8qKjI4LZWXThdvXoVlZWV8PHx0dru4+ODEydOGHSM33//HUuXLkW7du2k+VOLFy9G27ZtddomJSUhMTFRelxYWAh/f394eXnB1dW15k/kPtRqNRQKBby8vGT/ofP2BkKDBdLPKPDnWS88HytrOGZhTfmuD5hvy2POLYv5tjxz5dzBwcHgtlZdOJlC586dDe6CU6lUUKl0LylTKpVm+6FQKBRmPb4xYqKB9DPAlv8p8ff7j57WataU7/qA+bY85tyymG/LM0fOjTmWVb/Snp6esLGxQU5Ojtb2nJwc+Pr6mu28ycnJCA0NRWRkpNnOYY1i7kz/2rkPKCuXNxYiIiJrZNWFk729PSIiIrBlyxZpm1qtxpYtWxAdHW228yYkJCA9PR379u0z2zmsUdtHAC8PoOQWr64jIiLSR/bCqbi4GGlpadKVcRkZGUhLS8PFixcBAImJiZg7dy6+++47HD9+HG+//TZKSkqkq+zIdJTKquE6AFi/XdZQiIiIrJLsc5z279+P7t27S481k7Pj4+OxcOFCDBgwAHl5eRg/fjyuXLmC8PBwbNy4UWfCOJlGn6eAH9cD63cAH40CHOzljoiIiMh6yF44devWDUJUf+n7yJEjMXLkSAtFVDXHKTk5GZWVlRY7p7WIDgcaewHZecCWP4De3eSOiIiIyHrIPlRnjerrHCegarjuuV5VX3+/Rt5YiIiIrA0LJ9IxqE9VAfXHIeDUebmjISIish4snEhHEx+g5xNVX3/zk7yxEBERWRMWTqTXiIFVn39OBS5mVd+WiIiovmDhpEd9XQDzbuGtgK6RQKUa+Gyu3NEQERFZBxZOetTnyeF3+9dbVXOd1m0HfjsgdzRERETyY+FE9xXSHHi1b9XXYz8HbhTKGw8REZHcWDhRtf75BhDgB2TmAKM+Bsor5I6IiIhIPiycqFpOjsDsiYDKHtixFxj9CXC7TO6oiIiI5MHCSQ9ODtfWOhiY9RFgawOs3Qa89C5wJU/uqIiIiCyPhZMenByuq0c08N1ngJsLcOg40GMw8PX3QMktuSMjIiKyHBZOZLDOEcCaWUBYCFB8E/j8W6DDc8D7nwHrtwN5+XJHSEREZF6y3+SXapegJsCqZGDNVmD6d0DGZSBlY9UHADR0BZr7V00o93ADGrpVbWvgCKjsAHu7qvlS9naAnR2gUABKRdVnhQJQAIBCz/Z79tWEUAPX8m1QeAtQ8F8Gs6sr+a7p+00OQg1cu2aDotu1O+e1BfNteUIALip5Y1AIIYS8IVivwsJCuLm5oaCgAK6uriY/vlqtRm5uLry9vaFU1r6fOiGA/X9W9TbtTgNOnJM7IiIiqsvs7QR2Lcox+d9NY/7es8eJakyhACLbVn0AwM1bQEYmcPYikJ0HXC8ArhdWfb51GygtA8rKgdLyqq/LK6qKL6EGBO58LQC1+OtriKp9avWdbQ8TsBBQCwGl4iG6rchwdSDfte/fSgGhFlAoFbjTR0tmxXxbmp0VVC1WEIL1SU5ORnJyMiorK+UOpVZp4Fh1BV7rYLkj0U+tFnf18PGXnLkx35bHnFsW8215VTmXN4baNz5kAbyqjoiIiPRh4URERERkIBZORERERAZi4URERERkIBZORERERAZi4URERERkIBZORERERAZi4aRHcnIyQkNDERkZKXcoREREZEVYOOnBdZyIiIhIHxZORERERAZi4URERERkIN6rrhrizh0+CwsLzXJ8tVqNoqIiODg4mPQuz6Qf821ZzLflMeeWxXxbnrlyrvk7Lwy4szcLp2oUFRUBAPz9/WWOhIiIiMytqKgIbm5u1bZRCEPKq3pKrVYjKysLLi4uUChMf+frwsJC+Pv749KlS3B1dTX58Ukb821ZzLflMeeWxXxbnrlyLoRAUVER/Pz8HtiTxR6naiiVSjRt2tTs53F1deUPnQUx35bFfFsec25ZzLflmSPnD+pp0uCgLBEREZGBWDgRERERGYiFk4xUKhUmTJgAlUoldyj1AvNtWcy35THnlsV8W5415JyTw4mIiIgMxB4nIiIiIgOxcCIiIiIyEAsnIiIiIgOxcJJJcnIygoKC4ODggKioKOzdu1fukGqlnTt3ok+fPvDz84NCocCqVau09gshMH78eDRu3BiOjo6IiYnB6dOntdrk5+dj0KBBcHV1hbu7O4YNG4bi4mILPovaY8qUKYiMjISLiwu8vb3Rr18/nDx5UqvN7du3kZCQgEaNGsHZ2RnPP/88cnJytNpcvHgRvXv3RoMGDeDt7Y0xY8agoqLCkk+l1pg1axbatWsnrVsTHR2NDRs2SPuZb/P69NNPoVAo8O6770rbmHPT+uijj6BQKLQ+QkJCpP3Wlm8WTjJYunQpEhMTMWHCBBw8eBBhYWGIjY1Fbm6u3KHVOiUlJQgLC0NycrLe/VOnTsXMmTMxe/Zs7NmzB05OToiNjcXt27elNoMGDcKxY8eQmpqKdevWYefOnRg+fLilnkKtsmPHDiQkJOB///sfUlNTUV5ejl69eqGkpERqM3r0aKxduxYpKSnYsWMHsrKy8Nxzz0n7Kysr0bt3b5SVleGPP/7Ad999h4ULF2L8+PFyPCWr17RpU3z66ac4cOAA9u/fj6eeegp9+/bFsWPHADDf5rRv3z588803aNeundZ25tz0WrdujezsbOnjt99+k/ZZXb4FWVzHjh1FQkKC9LiyslL4+fmJKVOmyBhV7QdA/Pzzz9JjtVotfH19xeeffy5tu3HjhlCpVOLHH38UQgiRnp4uAIh9+/ZJbTZs2CAUCoXIzMy0WOy1VW5urgAgduzYIYSoyq+dnZ1ISUmR2hw/flwAELt37xZCCPHLL78IpVIprly5IrWZNWuWcHV1FaWlpZZ9ArVUw4YNxbx585hvMyoqKhItW7YUqampomvXruKdd94RQvA9bg4TJkwQYWFhevdZY77Z42RhZWVlOHDgAGJiYqRtSqUSMTEx2L17t4yR1T0ZGRm4cuWKVq7d3NwQFRUl5Xr37t1wd3dHhw4dpDYxMTFQKpXYs2ePxWOubQoKCgAAHh4eAIADBw6gvLxcK+chISEICAjQynnbtm3h4+MjtYmNjUVhYaHUi0L6VVZW4qeffkJJSQmio6OZbzNKSEhA7969tXIL8D1uLqdPn4afnx+aN2+OQYMG4eLFiwCsM9+8V52FXb16FZWVlVovMAD4+PjgxIkTMkVVN125cgUA9OZas+/KlSvw9vbW2m9rawsPDw+pDemnVqvx7rvv4oknnkCbNm0AVOXT3t4e7u7uWm3vzbm+10Szj3QdPXoU0dHRuH37NpydnfHzzz8jNDQUaWlpzLcZ/PTTTzh48CD27duns4/vcdOLiorCwoUL8eijjyI7OxsTJ05Ely5d8Oeff1plvlk4EVGNJCQk4M8//9Sai0Dm8eijjyItLQ0FBQVYvnw54uPjsWPHDrnDqpMuXbqEd955B6mpqXBwcJA7nHohLi5O+rpdu3aIiopCYGAgli1bBkdHRxkj049DdRbm6ekJGxsbnSsCcnJy4OvrK1NUdZMmn9Xl2tfXV2dSfkVFBfLz8/l6VGPkyJFYt24dtm3bhqZNm0rbfX19UVZWhhs3bmi1vzfn+l4TzT7SZW9vj+DgYERERGDKlCkICwvDjBkzmG8zOHDgAHJzc/HYY4/B1tYWtra22LFjB2bOnAlbW1v4+Pgw52bm7u6ORx55BGfOnLHK9zgLJwuzt7dHREQEtmzZIm1Tq9XYsmULoqOjZYys7mnWrBl8fX21cl1YWIg9e/ZIuY6OjsaNGzdw4MABqc3WrVuhVqsRFRVl8ZitnRACI0eOxM8//4ytW7eiWbNmWvsjIiJgZ2enlfOTJ0/i4sWLWjk/evSoVsGampoKV1dXhIaGWuaJ1HJqtRqlpaXMtxn06NEDR48eRVpamvTRoUMHDBo0SPqaOTev4uJinD17Fo0bN7bO97jJp5vTA/30009CpVKJhQsXivT0dDF8+HDh7u6udUUAGaaoqEgcOnRIHDp0SAAQ06ZNE4cOHRIXLlwQQgjx6aefCnd3d7F69Wpx5MgR0bdvX9GsWTNx69Yt6RhPP/20aN++vdizZ4/47bffRMuWLcXAgQPlekpW7e233xZubm5i+/btIjs7W/q4efOm1Oatt94SAQEBYuvWrWL//v0iOjpaREdHS/srKipEmzZtRK9evURaWprYuHGj8PLyEklJSXI8Jas3btw4sWPHDpGRkSGOHDkixo0bJxQKhdi0aZMQgvm2hLuvqhOCOTe19957T2zfvl1kZGSI33//XcTExAhPT0+Rm5srhLC+fLNwkslXX30lAgIChL29vejYsaP43//+J3dItdK2bdsEAJ2P+Ph4IUTVkgQffvih8PHxESqVSvTo0UOcPHlS6xjXrl0TAwcOFM7OzsLV1VUMGTJEFBUVyfBsrJ++XAMQCxYskNrcunVLjBgxQjRs2FA0aNBA9O/fX2RnZ2sd5/z58yIuLk44OjoKT09P8d5774ny8nILP5vaYejQoSIwMFDY29sLLy8v0aNHD6loEoL5toR7Cyfm3LQGDBggGjduLOzt7UWTJk3EgAEDxJkzZ6T91pZvhRBCmL4fi4iIiKju4RwnIiIiIgOxcCIiIiIyEAsnIiIiIgOxcCIiIiIyEAsnIiIiIgOxcCIiIiIyEAsnIiIiIgOxcCIiIiIyEAsnIiIiIgOxcCKiemnw4MHo16+f3GEQUS3DwomIiIjIQCyciKhOW758Odq2bQtHR0c0atQIMTExGDNmDL777jusXr0aCoUCCoUC27dvBwBcunQJL774Itzd3eHh4YG+ffvi/Pnz0vE0PVUTJ06El5cXXF1d8dZbb6GsrEyeJ0hEFmUrdwBEROaSnZ2NgQMHYurUqejfvz+Kioqwa9cuvPbaa7h48SIKCwuxYMECAICHhwfKy8sRGxuL6Oho7Nq1C7a2tpg8eTKefvppHDlyBPb29gCALVu2wMHBAdu3b8f58+cxZMgQNGrUCB9//LGcT5eILICFExHVWdnZ2aioqMBzzz2HwMBAAEDbtm0BAI6OjigtLYWvr6/U/vvvv4darca8efOgUCgAAAsWLIC7uzu2b9+OXr16AQDs7e0xf/58NGjQAK1bt8a///1vjBkzBpMmTYJSyY58orqMP+FEVGeFhYWhR48eaNu2Lf7+979j7ty5uH79+n3bHz58GGfOnIGLiwucnZ3h7OwMDw8P3L59G2fPntU6boMGDaTH0dHRKC4uxqVLl8z6fIhIfuxxIqI6y8bGBqmpqfjjjz+wadMmfPXVV/jggw+wZ88eve2Li4sRERGBH374QWefl5eXucMlolqAhRMR1WkKhQJPPPEEnnjiCYwfPx6BgYH4+eefYW9vj8rKSq22jz32GJYuXQpvb2+4urre95iHDx/GrVu34OjoCAD43//+B2dnZ/j7+5v1uRCR/DhUR0R11p49e/DJJ59g//79uHjxIlauXIm8vDy0atUKQUFBOHLkCE6ePImrV6+ivLwcgwYNgqenJ/r27Ytdu3YhIyMD27dvx6hRo3D58mXpuGVlZRg2bBjS09Pxyy+/YMKECRg5ciTnNxHVA+xxIqI6y9XVFTt37sT06dNRWFiIwMBAfPHFF4iLi0OHDh2wfft2dOjQAcXFxdi2bRu6deuGnTt34p///Ceee+45FBUVoUmTJujRo4dWD1SPHj3QsmVLPPnkkygtLcXAgQPx0UcfyfdEichiFEIIIXcQRES1xeDBg3Hjxg2sWrVK7lCISAbsVyYiIiIyEAsnIiIiIgNxqI6IiIjIQOxxIiIiIjIQCyciIiIiA7FwIiIiIjIQCyciIiIiA7FwIiIiIjIQCyciIiIiA7FwIiIiIjIQCyciIiIiA7FwIiIiIjLQ/wPrwS+xY6at3QAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# viz: loss curve on a log axis, the canonical training-curve view\n",
    "plt.figure(figsize=(6, 4))\n",
    "plt.plot(losses, color=\"#1E40FF\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(\"step\"); plt.ylabel(\"MSE loss (log scale)\")\n",
    "plt.title(f\"Gradient descent convergence ({STEPS} steps, lr=0.1)\")\n",
    "plt.grid(True, alpha=0.3)\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5420bf6",
   "metadata": {},
   "source": [
    "> **Interpretation.** The loss drops fast then flattens at the noise floor (it cannot go below the variance of the `0.1*noise` term). The recovered parameters sit within 0.1 of the truth. On a log axis the early steep descent and the late plateau are both visible at once, which is why loss curves are almost always plotted log-y.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef890fcf",
   "metadata": {},
   "source": [
    "### The same fit, two more ways: the implementation ladder\n",
    "\n",
    "The spec's implementation ladder says: derive the math, write it from scratch, call the library, and assert they agree. We just did from-scratch. Now `torch.autograd` with an explicit loop, and the closed-form normal equation. All three recover the same parameters from the same data, which is the agreement check that proves none of them lied.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "ae5732e0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:02.962578Z",
     "iopub.status.busy": "2026-06-10T18:44:02.962504Z",
     "iopub.status.idle": "2026-06-10T18:44:03.008006Z",
     "shell.execute_reply": "2026-06-10T18:44:03.007549Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch:  w=2.9931, b=4.9913\n",
      "normal: w=2.9931, b=4.9913\n",
      "[ ok ] from-scratch, torch, and normal equation all agree\n"
     ]
    }
   ],
   "source": [
    "# torch path: autograd does the gradient we derived by hand\n",
    "torch.manual_seed(SEED)\n",
    "Xt = torch.tensor(Xreg, dtype=torch.float32)\n",
    "yt = torch.tensor(yreg, dtype=torch.float32)\n",
    "wt = torch.zeros(1, requires_grad=True)\n",
    "bt = torch.zeros(1, requires_grad=True)\n",
    "for step in range(STEPS):\n",
    "    pred = (Xt @ wt + bt).squeeze()        # forward\n",
    "    loss = 0.5 * ((pred - yt) ** 2).mean()  # MSE\n",
    "    loss.backward()                         # backward\n",
    "    with torch.no_grad():                   # update without tracking\n",
    "        wt -= 0.1 * wt.grad\n",
    "        bt -= 0.1 * bt.grad\n",
    "        wt.grad = None; bt.grad = None      # zero grads the modern way (not .data)\n",
    "print(f\"torch:  w={wt.item():.4f}, b={bt.item():.4f}\")\n",
    "\n",
    "# closed form: normal equation w = (A^T A)^-1 A^T y on the bias-augmented design\n",
    "A = np.hstack([Xreg, np.ones((Xreg.shape[0], 1))])   # append a ones column for the bias\n",
    "theta = np.linalg.solve(A.T @ A, A.T @ yreg)         # solve, never invert explicitly\n",
    "print(f\"normal: w={theta[0]:.4f}, b={theta[1]:.4f}\")\n",
    "\n",
    "assert abs(wt.item() - w_hat[0]) < 0.05, \"torch GD should match our from-scratch GD\"\n",
    "assert abs(theta[0] - 3.0) < 0.05 and abs(theta[1] - 5.0) < 0.05, \"normal equation should recover (3, 5)\"\n",
    "print(\"[ ok ] from-scratch, torch, and normal equation all agree\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa3dba4e",
   "metadata": {},
   "source": [
    "> **Note:** we zero gradients with `wt.grad = None` inside `torch.no_grad()`, never by mutating `.data`. PyTorch accumulates gradients across `backward()` calls, so a loop that forgets to clear them sums every step's gradient, which is the gradient-accumulation footgun the deep-learning chapters demonstrate in full. We solve the normal equation with `np.linalg.solve`, not by forming `(A^T A)^{-1}` explicitly, because solving is both faster and more numerically stable than inverting.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bffa21e",
   "metadata": {},
   "source": [
    "### Experiment log\n",
    "\n",
    "The Karpathy convention: a table recording every named change with its measured numbers. It doubles as your expected-value reference. If your run does not produce roughly these numbers, something is wrong.\n",
    "\n",
    "| method | lr | steps (full / FAST) | final w | final b | final loss |\n",
    "|---|---|---|---|---|---|\n",
    "| from-scratch GD | 0.1 | 5000 / 500 | ~3.00 | ~5.00 | ~0.0049 |\n",
    "| torch autograd GD | 0.1 | 5000 / 500 | ~3.00 | ~5.00 | matches GD |\n",
    "| normal equation | closed form | n/a | ~3.00 | ~5.00 | global minimum |\n",
    "\n",
    "> **Key takeaways.** Gradient descent is four lines: forward, backward, update, track. Ground-truth synthesis lets you assert recovery, the strongest possible check. The from-scratch, autograd, and closed-form paths converge to the same answer, which is the agreement check that retires any doubt about the implementation. Loss curves go on a log y-axis.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb9f064f",
   "metadata": {},
   "source": [
    "## Safety lens\n",
    "\n",
    "The math is value-neutral; the practice has failure modes that matter for safety. Two are about numerical stability hiding bugs that look like model failures, one is about reproducibility as a safety property.\n",
    "\n",
    "You already met the first: a naive softmax returns `nan` on a large logit, and an evaluation pipeline that silently produces `nan` can log \"model failed safely\" when the truth is \"we never tested anything\". The same trap lives in variance estimated as `E[X^2] - E[X]^2` on near-equal scales, where catastrophic cancellation can return a negative variance. We demonstrate it, then show the stable two-pass estimator NumPy uses.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "0c800e37",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:03.008938Z",
     "iopub.status.busy": "2026-06-10T18:44:03.008847Z",
     "iopub.status.idle": "2026-06-10T18:44:03.011399Z",
     "shell.execute_reply": "2026-06-10T18:44:03.011045Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "naive  E[X^2]-E[X]^2 variance: 0.0000\n",
      "stable two-pass variance:      0.9813   (true spread is ~1.0)\n",
      "The naive formula loses almost all precision here; it can even go negative.\n"
     ]
    }
   ],
   "source": [
    "# deeper: catastrophic cancellation in the textbook variance formula\n",
    "big = 1e8 + rng.standard_normal(1000)        # values clustered around 1e8 with unit-scale spread\n",
    "naive_var = (big**2).mean() - big.mean()**2  # E[X^2] - E[X]^2 in float64\n",
    "stable_var = big.var()                        # NumPy's two-pass: mean((x - mean)^2)\n",
    "print(f\"naive  E[X^2]-E[X]^2 variance: {naive_var:.4f}\")\n",
    "print(f\"stable two-pass variance:      {stable_var:.4f}   (true spread is ~1.0)\")\n",
    "print(\"The naive formula loses almost all precision here; it can even go negative.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e93346ba",
   "metadata": {},
   "source": [
    "> **Caveat:** depending on your BLAS, `naive_var` above may print a wildly wrong or even negative number while `stable_var` stays near 1.0. That is the whole point: the algebra is identical, the floating-point behavior is not. Reach for `np.var` / `torch.var`, not the hand formula, exactly as you reach for `F.cross_entropy` over softmax-then-log.\n",
    "\n",
    "Reproducibility is the second safety property. Alignment research makes empirical claims like \"method X reduces harmful behavior by Y%\", and an unreproducible claim is an unfalsifiable one. The mechanical fixes are the ones this notebook already practices: seed every RNG (`np.random.default_rng`, `torch.manual_seed`, `random.seed`), re-seed before stochastic cells so isolated re-runs reproduce, and quote numbers only from seeded runs. The deeper habit underneath every safety practice in this curriculum: believe nothing until it has been measured. Every `assert` and every `check` in this notebook is that habit made executable.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1ce98ba",
   "metadata": {},
   "source": [
    "## Test yourself\n",
    "\n",
    "Three parts: concept self-checks with folded answers, two auto-checked problems, and a capstone with a rubric and a folded reference. Try each before you peek. Every answer is in this notebook; if unsure, re-run that section.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ef3d070",
   "metadata": {},
   "source": [
    "### Part A — Concepts\n",
    "\n",
    "1. You have `X` of shape `(3, 4)`. What shape is `X.sum(axis=0)`, and what is each entry? <details><summary>Answer</summary>Shape `(4,)`. Axis 0 (the 3 rows) is collapsed, leaving one sum per column. Re-run the reductions cell in Part 1 to see it.</details>\n",
    "2. Why does `X.mean(axis=1)` fail to broadcast back against `X` for row-centering, while `X.mean(axis=1, keepdims=True)` works? <details><summary>Answer</summary>Without `keepdims` the mean has shape `(n,)`, which aligns from the right against `X`'s columns, not its rows. `keepdims=True` gives `(n, 1)`, which stretches across the columns and subtracts from each row. This is Exercise 0.1.</details>\n",
    "3. In `(3, 4) @ (4, 5)`, which numbers must match and what is the output shape? <details><summary>Answer</summary>The inner dimensions (the two 4s) must match; they vanish. The output is the outer dimensions, `(3, 5)`.</details>\n",
    "4. The finite-difference cell in Part 3 showed the estimate getting better as `h` shrank, then it would get worse again past `1e-8` or so. Why both directions? <details><summary>Answer</summary>Larger `h` means the slope is measured over a curved interval (truncation error, falls as you shrink `h`). Too-small `h` means `f(x+h) - f(x-h)` subtracts two nearly equal floats and loses precision (round-off error, grows as you shrink `h`). The sweet spot is around `1e-5` for a centered difference.</details>\n",
    "5. A 99%-sensitive test for a 1%-prevalence disease comes back positive. Roughly what is the probability you have the disease, and why so low? <details><summary>Answer</summary>About 17% (Exercise 0.4). With 99 healthy people per sick one, the 5% false-positive rate produces more positives than the disease itself. The base rate dominates.</details>\n",
    "6. Why do we minimize negative log-likelihood instead of maximizing the likelihood directly? <details><summary>Answer</summary>The likelihood is a product of many probabilities, which underflows to exactly 0 in floating point (the 2000-event cell in Part 4). The log turns the product into a sum that stays finite, and maximizing the log is the same as maximizing the original. The minus is just so we can minimize.</details>\n",
    "7. Why subtract the max logit before the softmax exponentiation, and why is it exact? <details><summary>Answer</summary>`exp` of a large logit overflows to `inf` and `inf/inf` is `nan` (the deliberate-failure cell). Subtracting a constant from every logit multiplies numerator and denominator by the same factor, which cancels, so the softmax output is unchanged. It is the log-sum-exp trick.</details>\n",
    "8. Is KL divergence a distance? <details><summary>Answer</summary>No. It is non-negative and zero only when the distributions are equal, but it is asymmetric (`KL(p||q) != KL(q||p)`) and violates the triangle inequality, so it is not a metric. Part 5 checks the asymmetry indirectly via the identity `KL = H(p,q) - H(p)`.</details>\n",
    "9. The Part 6 training loop runs a fixed step count with no early `break`. Why is that safer than breaking when the loss \"stops changing\"? <details><summary>Answer</summary>An early-break sentinel can fire on the first iteration before any update happens, if the stopping condition happens to hold at initialization. The named k-means early-break bug is exactly this. A fixed step count cannot exit spuriously.</details>\n",
    "10. You hand-derive a gradient and want to trust it. What is the cheapest check, and what is the strongest? <details><summary>Answer</summary>Cheapest: a centered finite-difference gradient check (Exercise 0.3), which needs only forward evaluations. Strongest: `torch.autograd` on the same expression with `torch.testing.assert_close` (Part 3), which checks the exact analytic gradient, not an approximation.</details>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67b7bf9e",
   "metadata": {},
   "source": [
    "### Part B — Auto-checked problems\n",
    "\n",
    "**Problem B1 — Cosine similarity** · `Difficulty 2/5 · ~8 min`\n",
    "\n",
    "Implement `cosine_similarity(x, y)` = $\\frac{\\mathbf{x}^\\top\\mathbf{y}}{\\|\\mathbf{x}\\|\\,\\|\\mathbf{y}\\|}$ for two 1D vectors. The check asserts a vector is perfectly similar to itself (1.0), that orthogonal vectors give 0, and that opposite vectors give -1.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "90576bfc",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:03.012511Z",
     "iopub.status.busy": "2026-06-10T18:44:03.012401Z",
     "iopub.status.idle": "2026-06-10T18:44:03.017256Z",
     "shell.execute_reply": "2026-06-10T18:44:03.016857Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ -- ] B1 cosine_similarity: not attempted yet — fill in the TODO above, then re-run.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def cosine_similarity(x, y):\n",
    "    x = np.asarray(x, dtype=float); y = np.asarray(y, dtype=float)\n",
    "    # TODO: dot(x, y) divided by the product of the two L2 norms\n",
    "    result = None\n",
    "    attempted(result)\n",
    "    return result\n",
    "\n",
    "def _cosine():\n",
    "    v = np.array([1.0, 2.0, 3.0])\n",
    "    check_close(cosine_similarity(v, v), 1.0, msg=\"a vector is perfectly aligned with itself\")\n",
    "    check_close(cosine_similarity([1.0, 0.0], [0.0, 1.0]), 0.0, msg=\"orthogonal vectors -> 0\")\n",
    "    check_close(cosine_similarity([1.0, 0.0], [-1.0, 0.0]), -1.0, msg=\"opposite vectors -> -1\")\n",
    "\n",
    "check(\"B1 cosine_similarity\", _cosine)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e230219",
   "metadata": {},
   "source": [
    "<details><summary>Hint 1 (conceptual)</summary>The numerator is `x @ y`. The denominator is `np.linalg.norm(x) * np.linalg.norm(y)`.</details>\n",
    "\n",
    "<details><summary>Hint 2 (the line)</summary>`result = (x @ y) / (np.linalg.norm(x) * np.linalg.norm(y))`</details>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "b2be871f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:03.018391Z",
     "iopub.status.busy": "2026-06-10T18:44:03.018303Z",
     "iopub.status.idle": "2026-06-10T18:44:03.021664Z",
     "shell.execute_reply": "2026-06-10T18:44:03.021365Z"
    },
    "collapsed": true,
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ ok ] B1 cosine_similarity\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#@title Solution { display-mode: \"form\" }\n",
    "# solution: redefines cosine_similarity; the check below re-verifies the reference.\n",
    "def cosine_similarity(x, y):\n",
    "    x = np.asarray(x, dtype=float); y = np.asarray(y, dtype=float)\n",
    "    return (x @ y) / (np.linalg.norm(x) * np.linalg.norm(y))\n",
    "\n",
    "check(\"B1 cosine_similarity\", _cosine, required=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66bc21a6",
   "metadata": {},
   "source": [
    "**Problem B2 — Stable log-sum-exp** · `Difficulty 3/5 · ~12 min`\n",
    "\n",
    "The denominator of a log-softmax is $\\log\\sum_k e^{z_k}$, which overflows for the same reason softmax does. Implement `logsumexp(z)` stably by factoring out the max: $\\log\\sum_k e^{z_k} = m + \\log\\sum_k e^{z_k - m}$ where $m = \\max_k z_k$. The check confirms it matches the naive formula on small inputs and stays finite on a logit of 1000.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "2ad95b47",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:03.022526Z",
     "iopub.status.busy": "2026-06-10T18:44:03.022445Z",
     "iopub.status.idle": "2026-06-10T18:44:03.025903Z",
     "shell.execute_reply": "2026-06-10T18:44:03.025641Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ -- ] B2 logsumexp: not attempted yet — fill in the TODO above, then re-run.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def logsumexp(z):\n",
    "    z = np.asarray(z, dtype=float)\n",
    "    # TODO 1: m = the max of z\n",
    "    m = None\n",
    "    attempted(m)\n",
    "    # TODO 2: return m + log(sum(exp(z - m)))\n",
    "    return m + np.log(np.exp(z - m).sum())\n",
    "\n",
    "def _lse():\n",
    "    small = np.array([0.1, 0.2, 0.3])\n",
    "    naive = np.log(np.exp(small).sum())                     # safe for small values\n",
    "    check_close(logsumexp(small), naive, msg=\"logsumexp must match the naive value on small inputs\")\n",
    "    big = logsumexp(np.array([1.0, 2.0, 1000.0]))\n",
    "    assert np.isfinite(big) and abs(big - 1000.0) < 1e-6, \\\n",
    "        \"logsumexp of a 1000 logit should be ~1000 and finite, not inf\"\n",
    "\n",
    "check(\"B2 logsumexp\", _lse)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05671cd1",
   "metadata": {},
   "source": [
    "<details><summary>Hint 1 (conceptual)</summary>Subtract the max `m` before exponentiating, then add `m` back outside the log. This is the same trick as the stable softmax, written for the log of the sum.</details>\n",
    "\n",
    "<details><summary>Help — \"logsumexp returns inf\"</summary>You exponentiated `z` directly instead of `z - m`. `exp(1000)` is `inf`. Subtract the max first.</details>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "7ed8ee38",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-06-10T18:44:03.027165Z",
     "iopub.status.busy": "2026-06-10T18:44:03.027082Z",
     "iopub.status.idle": "2026-06-10T18:44:03.029977Z",
     "shell.execute_reply": "2026-06-10T18:44:03.029577Z"
    },
    "collapsed": true,
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ ok ] B2 logsumexp\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#@title Solution { display-mode: \"form\" }\n",
    "# solution: redefines logsumexp; the check below re-verifies the reference.\n",
    "def logsumexp(z):\n",
    "    z = np.asarray(z, dtype=float)\n",
    "    m = z.max()\n",
    "    return float(m + np.log(np.exp(z - m).sum()))\n",
    "\n",
    "check(\"B2 logsumexp\", _lse, required=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df4e9514",
   "metadata": {},
   "source": [
    "### Part C — Capstone: gradient descent three ways, agreement-checked\n",
    "\n",
    "Redo the Part 6 fit as a self-contained project on a fresh two-feature problem `y = 2*x1 - 3*x2 + 1 + noise`, and prove three implementations agree.\n",
    "\n",
    "**Deliverables**\n",
    "1. Synthesize 200 points from the law above with a seeded `default_rng`, true parameters `w = [2, -3]`, `b = 1`.\n",
    "2. Fit it three ways: from-scratch NumPy gradient descent (reuse `mse_loss` / `mse_gradient` generalized to two features), the closed-form normal equation, and `sklearn.linear_model.LinearRegression`.\n",
    "3. Overlay the loss curve of the from-scratch run on a log y-axis.\n",
    "\n",
    "**Self-assessment (pass / partial / fail)**\n",
    "- (a) all three recover `(2, -3, 1)` to two decimals;\n",
    "- (b) the from-scratch and normal-equation weights agree to within 0.05;\n",
    "- (c) the loss curve is plotted with a log y-axis and labeled axes;\n",
    "- (d) you re-seeded the RNG so the run reproduces;\n",
    "- (e) you can state in one sentence why all three converge to the same weights (the loss is convex, so there is one global minimum).\n",
    "\n",
    "<details><summary>My solution (reference, runs in well under a second on CPU)</summary>\n",
    "\n",
    "```python\n",
    "gen = np.random.default_rng(SEED)\n",
    "n = 200\n",
    "Xc = gen.standard_normal((n, 2))\n",
    "true_w = np.array([2.0, -3.0]); true_b = 1.0\n",
    "yc = Xc @ true_w + true_b + 0.1 * gen.standard_normal(n)\n",
    "\n",
    "# from-scratch GD (mse_loss / mse_gradient already handle vector w)\n",
    "w = np.zeros(2); b = 0.0; hist = []\n",
    "for _ in range(STEPS):\n",
    "    hist.append(mse_loss(w, b, Xc, yc))\n",
    "    dw, db = mse_gradient(w, b, Xc, yc)\n",
    "    w = w - 0.1 * dw; b = b - 0.1 * db\n",
    "print(\"GD     \", w.round(3), round(b, 3))\n",
    "\n",
    "# normal equation on the bias-augmented design\n",
    "A = np.hstack([Xc, np.ones((n, 1))])\n",
    "theta = np.linalg.solve(A.T @ A, A.T @ yc)\n",
    "print(\"normal \", theta[:2].round(3), round(theta[2], 3))\n",
    "\n",
    "# sklearn\n",
    "from sklearn.linear_model import LinearRegression\n",
    "m = LinearRegression().fit(Xc, yc)\n",
    "print(\"sklearn\", m.coef_.round(3), round(float(m.intercept_), 3))\n",
    "\n",
    "assert np.allclose(w, true_w, atol=0.05) and abs(b - true_b) < 0.05\n",
    "assert np.allclose(theta[:2], true_w, atol=0.05)\n",
    "assert np.allclose(m.coef_, true_w, atol=0.05)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(hist); plt.yscale(\"log\"); plt.xlabel(\"step\"); plt.ylabel(\"MSE (log)\")\n",
    "plt.title(\"capstone: GD convergence\"); plt.grid(True, alpha=0.3); plt.show()\n",
    "```\n",
    "All three land on `(2, -3, 1)` because squared-error loss on a full-rank design is convex with a single global minimum; GD descends to it, the normal equation solves for it directly, and sklearn solves the same equation.</details>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7342c47",
   "metadata": {},
   "source": [
    "## Reflection\n",
    "\n",
    "Write ~150 words, in the cell below, on the dumbest bug you hit in this notebook and how you found it. Maybe it was a `keepdims` you forgot and a broadcast that silently did the wrong thing, or a softmax that returned `nan`, or a gradient your finite-difference checker flagged. Nobody grades this. Writing it is the point: the skill that carries you through the next 26 chapters is not the math, it is the habit of noticing the gap between what you expected and what the cell printed, and chasing it down. The whole curriculum is built on that loop. State plainly which check first told you something was wrong, and what you printed to confirm the fix.\n",
    "\n",
    "*(Replace this text with your reflection.)*\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67a02c8b",
   "metadata": {},
   "source": [
    "## Going further\n",
    "\n",
    "- 3Blue1Brown, *Essence of Linear Algebra* and *Essence of Calculus* (YouTube): the geometric intuition behind every operation in Parts 2 and 3.\n",
    "- Karpathy, *The spelled-out intro to neural networks and backpropagation: building micrograd*: the finite-difference opener of Part 3 expanded into a full autograd engine.\n",
    "- Géron, *Hands-On ML* 3e, the standalone `math_linear_algebra` and `math_differential_calculus` notebooks: worked examples for everything here.\n",
    "- d2l.ai, *Appendix: Mathematics for Deep Learning*: the densest reference for the DL-specific math, especially the information-theory and maximum-likelihood sections.\n",
    "- MacKay, *Information Theory, Inference, and Learning Algorithms*, ch. 1-2: the entropy intuitions of Part 5, less computational and more philosophical.\n",
    "- The CS229 linear algebra and probability cheat sheets (Amidi): two pages each, worth printing.\n",
    "\n",
    "## What this enables\n",
    "\n",
    "- **Ch 01 — The ML Landscape**: every algorithm there is now \"define a loss, take its gradient, take a step\", the four bullets you implemented in Part 6.\n",
    "- **Ch 03 — Classification**: cross-entropy, softmax, and the base-rate fallacy from Parts 4 and 5 are that chapter's entire toolkit.\n",
    "- **Ch 09-15 — the deep-learning arc**: the `X @ W.T + b` layer, the grad-check, the stable softmax, and the four-line training loop are the atoms every transformer is built from.\n",
    "\n",
    "> **The gap this leaves.** We fit one convex problem with a hand-written loop. The next chapters stack thousands of these layers into a function with millions of parameters and a loss landscape that is not convex, where the same gradient descent still works but initialization, learning-rate schedules, and architecture start to matter enormously. The math does not change. The engineering around it is the rest of the curriculum.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24acc3bc",
   "metadata": {},
   "source": [
    "---\n",
    "*Built top-to-bottom. If every check above printed `[ ok ]`, you have reproduced the chapter. Total running time and last-verified date written by CI.*\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "obvix-nb",
   "language": "python",
   "name": "obvix-nb"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  },
  "obvix": {
   "title": "Ch 00 — Math & Python Prereqs"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}