{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "47a2dcb6",
   "metadata": {},
   "source": [
    "# Lesson 30: Neural Network Fundamentals — Learning the Projection\n",
    "\n",
    "Lesson 29 hand-picked a projection direction $w$ (the difference of class means) to separate two classes. That worked, but it required a human to notice a good heuristic. This lesson replaces the human with **gradient descent**: an automatic procedure that *learns* $w$ and $b$ directly from data by repeatedly nudging them to reduce a loss. The recipe — forward pass, loss, backward pass, update — is the entire training loop behind every neural network in this course, no matter how large."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19a1b6b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62b9aea2",
   "metadata": {},
   "source": [
    "## The same dataset as Lesson 29"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb7a17b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(1)\n",
    "class_a = rng.normal(loc=[-2, -1], scale=0.8, size=(60, 2))\n",
    "class_b = rng.normal(loc=[2, 1.5], scale=0.8, size=(60, 2))\n",
    "X = np.vstack([class_a, class_b])\n",
    "y = np.concatenate([np.zeros(60), np.ones(60)])  # class A = 0, class B = 1\n",
    "\n",
    "plt.scatter(*class_a.T, s=15, label='class A (y=0)')\n",
    "plt.scatter(*class_b.T, s=15, label='class B (y=1)')\n",
    "plt.legend(fontsize=8)\n",
    "plt.gca().set_aspect('equal')\n",
    "plt.title('Same two-class dataset as Lesson 29')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e100244f",
   "metadata": {},
   "source": [
    "## From a raw score to a trainable loss\n",
    "\n",
    "Classification accuracy (right or wrong) is flat almost everywhere and jumps discontinuously at the decision boundary — it has no useful gradient to follow. Instead, squash the raw score $z=w^\\top x+b$ through the **sigmoid** function to get a probability, and measure error with **binary cross-entropy**:\n",
    "\n",
    "$$\\sigma(z) = \\frac{1}{1+e^{-z}}, \\qquad L = -\\frac{1}{n}\\sum_i \\big[y_i \\log \\sigma(z_i) + (1-y_i)\\log(1-\\sigma(z_i))\\big]$$\n",
    "\n",
    "Both pieces are smooth, so $L$ has a well-defined gradient everywhere — the loss decreases smoothly as predictions get closer to being correct, instead of only changing at the boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a497ef8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(z):\n",
    "    return 1 / (1 + np.exp(-z))\n",
    "\n",
    "zs = np.linspace(-6, 6, 200)\n",
    "plt.plot(zs, sigmoid(zs))\n",
    "plt.axhline(0.5, color='gray', linestyle='--', linewidth=1)\n",
    "plt.axvline(0, color='gray', linestyle='--', linewidth=1)\n",
    "plt.xlabel('z = w.x + b')\n",
    "plt.ylabel('sigma(z)')\n",
    "plt.title('The sigmoid activation')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc6f5345",
   "metadata": {},
   "source": [
    "## Backpropagation: the chain rule, applied\n",
    "\n",
    "\"Backprop\" is just the chain rule from Lessons 12-13, applied to a composed function: $L$ depends on $\\sigma(z)$, which depends on $z=w^\\top x + b$, which depends on $w$ and $b$. Working backward through that chain (skipping the algebra, which is a standard but slightly tedious simplification) gives a remarkably clean result:\n",
    "\n",
    "$$\\frac{\\partial L}{\\partial z_i} = \\sigma(z_i) - y_i, \\qquad \\frac{\\partial L}{\\partial w} = \\frac{1}{n}X^\\top(\\sigma(z)-y), \\qquad \\frac{\\partial L}{\\partial b} = \\frac{1}{n}\\sum_i(\\sigma(z_i)-y_i)$$\n",
    "\n",
    "In words: the gradient is just the *prediction error*, averaged (weighted by $x$ for $w$'s gradient). Wildly wrong predictions push the weights hard; correct ones barely move them at all."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86fec90d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def forward(X, w, b):\n",
    "    z = X @ w + b\n",
    "    return sigmoid(z), z\n",
    "\n",
    "def bce_loss(p, y, eps=1e-9):\n",
    "    return -np.mean(y * np.log(p + eps) + (1 - y) * np.log(1 - p + eps))\n",
    "\n",
    "def backward(X, p, y):\n",
    "    grad_z = (p - y) / len(y)\n",
    "    grad_w = X.T @ grad_z\n",
    "    grad_b = grad_z.sum()\n",
    "    return grad_w, grad_b"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7328979",
   "metadata": {},
   "source": [
    "### Sanity check: gradient checking\n",
    "\n",
    "Before trusting any hand-derived gradient, it's standard practice to check it numerically: perturb each parameter by a tiny amount and see how much the loss actually changes, then compare to the analytical formula."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1d1d4f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "w_test, b_test = np.array([0.5, -0.3]), 0.2\n",
    "p_test, _ = forward(X, w_test, b_test)\n",
    "grad_w_analytic, grad_b_analytic = backward(X, p_test, y)\n",
    "\n",
    "def loss_at(w, b):\n",
    "    p, _ = forward(X, w, b)\n",
    "    return bce_loss(p, y)\n",
    "\n",
    "eps = 1e-5\n",
    "grad_w_numeric = np.zeros(2)\n",
    "for i in range(2):\n",
    "    w_plus, w_minus = w_test.copy(), w_test.copy()\n",
    "    w_plus[i] += eps\n",
    "    w_minus[i] -= eps\n",
    "    grad_w_numeric[i] = (loss_at(w_plus, b_test) - loss_at(w_minus, b_test)) / (2 * eps)\n",
    "grad_b_numeric = (loss_at(w_test, b_test + eps) - loss_at(w_test, b_test - eps)) / (2 * eps)\n",
    "\n",
    "print(f'analytic grad_w: {grad_w_analytic}, numeric: {grad_w_numeric}')\n",
    "print(f'analytic grad_b: {grad_b_analytic:.6f}, numeric: {grad_b_numeric:.6f}')\n",
    "print(f'max discrepancy: {max(np.abs(grad_w_analytic - grad_w_numeric).max(), abs(grad_b_analytic - grad_b_numeric)):.2e}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26944bc6",
   "source": "### Aside: what's a tensor?\n\nThis is the first appearance of PyTorch in the course, and with it, the word **tensor**. A tensor is just the general term for a grid of numbers of any dimensionality: a scalar is a 0-dimensional tensor, a vector is 1-dimensional, a matrix is 2-dimensional, and a stack of matrices (e.g. a batch of RGB images, indexed by `[batch, channel, height, width]`) is a 4-dimensional tensor. Every array used so far in this course — NumPy arrays, images, weight matrices — has really been a tensor all along; PyTorch's `torch.Tensor` is simply NumPy's `ndarray` with two extras bolted on: it can track a `.grad` for automatic differentiation (used below), and it can live on a GPU instead of the CPU for fast, parallel computation. A tensor's `.shape` (e.g. `(100, 2)` for `X_t` below) says exactly what NumPy's `.shape` would say for the same array — the two libraries are close enough that converting between them (`torch.tensor(numpy_array)`, `tensor.numpy()`) is essentially free.",
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "id": "cf7afe60",
   "metadata": {},
   "source": [
    "### Sanity check: PyTorch autograd\n",
    "\n",
    "As a second, independent check, we let PyTorch compute the same gradient automatically via `.backward()`, using its built-in `binary_cross_entropy_with_logits` (which combines the sigmoid and the loss in one numerically stable function)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b6a38e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_t = torch.tensor(X, dtype=torch.float64)\n",
    "y_t = torch.tensor(y, dtype=torch.float64)\n",
    "w_t = torch.tensor(w_test, dtype=torch.float64, requires_grad=True)\n",
    "b_t = torch.tensor(b_test, dtype=torch.float64, requires_grad=True)\n",
    "\n",
    "z_t = X_t @ w_t + b_t\n",
    "loss_t = torch.nn.functional.binary_cross_entropy_with_logits(z_t, y_t)\n",
    "loss_t.backward()\n",
    "\n",
    "print(f'our analytic grad_w: {grad_w_analytic}')\n",
    "print(f'torch autograd grad_w: {w_t.grad.numpy()}')\n",
    "print(f'max discrepancy: {np.abs(grad_w_analytic - w_t.grad.numpy()).max():.2e}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90ee69aa",
   "metadata": {},
   "source": [
    "Two independent checks — numerical finite differences and PyTorch's automatic differentiation — both agree with the hand-derived formula to many decimal places. This kind of double-checking is standard practice whenever you derive a gradient by hand."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "036dcdaa",
   "metadata": {},
   "source": [
    "## Training: the full loop\n",
    "\n",
    "Initialize $w, b$ randomly (not with a clever heuristic this time), then repeat: forward pass, compute loss, backward pass, take a small step *against* the gradient (gradient *descent*, since the gradient points toward increasing loss)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5353354e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def train(X, y, n_epochs=500, lr=0.5, seed=0):\n",
    "    rng_local = np.random.default_rng(seed)\n",
    "    w = rng_local.normal(size=X.shape[1]) * 0.1\n",
    "    b = 0.0\n",
    "    losses = []\n",
    "    for _ in range(n_epochs):\n",
    "        p, _ = forward(X, w, b)\n",
    "        losses.append(bce_loss(p, y))\n",
    "        grad_w, grad_b = backward(X, p, y)\n",
    "        w -= lr * grad_w\n",
    "        b -= lr * grad_b\n",
    "    return w, b, losses\n",
    "\n",
    "w_learned, b_learned, losses = train(X, y)\n",
    "final_pred = (sigmoid(X @ w_learned + b_learned) > 0.5).astype(float)\n",
    "print(f'learned w = {np.round(w_learned, 3)}, b = {b_learned:.3f}')\n",
    "print(f'final accuracy: {(final_pred == y).mean():.1%}')\n",
    "\n",
    "plt.plot(losses)\n",
    "plt.xlabel('epoch')\n",
    "plt.ylabel('loss')\n",
    "plt.title('Training loss')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3f7187e",
   "metadata": {},
   "source": [
    "### Compare to Lesson 29's hand-picked direction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf9c755a",
   "metadata": {},
   "outputs": [],
   "source": [
    "w_handpicked = class_b.mean(axis=0) - class_a.mean(axis=0)\n",
    "w_handpicked /= np.linalg.norm(w_handpicked)\n",
    "w_learned_normalized = w_learned / np.linalg.norm(w_learned)\n",
    "\n",
    "cos_angle = w_handpicked @ w_learned_normalized\n",
    "print(f'hand-picked direction (Lesson 29): {np.round(w_handpicked, 3)}')\n",
    "print(f'learned direction (this lesson):   {np.round(w_learned_normalized, 3)}')\n",
    "print(f'angle between them: {np.degrees(np.arccos(np.clip(cos_angle, -1, 1))):.1f} degrees')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6dc937b",
   "metadata": {},
   "source": [
    "Gradient descent, starting from nothing but random noise, rediscovers essentially the same direction a human picked by reasoning about class means — reassuring, but also a preview of the limits of a single neuron: it can only ever *rediscover* what a single linear projection is capable of."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e456fe81",
   "metadata": {},
   "source": [
    "## Training on the unsolvable problem\n",
    "\n",
    "Lesson 29 showed that *no* linear projection separates a ring from the disk it surrounds — the best exhaustive search over directions and thresholds found was 74% accuracy. What happens when gradient descent, rather than brute-force search, tries to solve the same problem?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c0de831",
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_in = rng.uniform(0, 2 * np.pi, 60)\n",
    "inner = np.stack([0.5 * np.cos(theta_in), 0.5 * np.sin(theta_in)], axis=1) + rng.normal(0, 0.1, (60, 2))\n",
    "theta_out = rng.uniform(0, 2 * np.pi, 60)\n",
    "outer = np.stack([2.0 * np.cos(theta_out), 2.0 * np.sin(theta_out)], axis=1) + rng.normal(0, 0.15, (60, 2))\n",
    "X_ring = np.vstack([inner, outer])\n",
    "y_ring = np.concatenate([np.zeros(60), np.ones(60)])\n",
    "\n",
    "w_ring, b_ring, losses_ring = train(X_ring, y_ring, n_epochs=2000)\n",
    "pred_ring = (sigmoid(X_ring @ w_ring + b_ring) > 0.5).astype(float)\n",
    "print(f'learned w = {np.round(w_ring, 4)} (norm={np.linalg.norm(w_ring):.4f}), b = {b_ring:.4f}')\n",
    "print(f'trained single-neuron accuracy: {(pred_ring == y_ring).mean():.1%}')\n",
    "print(f'(Lesson 29\\'s exhaustive-search ceiling for ANY linear separator was 74%)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d99055a",
   "metadata": {},
   "source": [
    "Gradient descent actually does *worse* here than brute-force search — it converges to a near-zero weight vector and roughly chance-level accuracy, instead of finding the small, lopsided arc that let a brute-force search eke out 74%. The dataset is (approximately) symmetric around the origin, so the *average* gradient pull from all the training points nearly cancels out in every direction, and the optimizer settles near $w\\approx 0$ rather than hunting for an asymmetric corner-case solution. Different failure mode, same underlying truth: a single linear projection cannot solve this problem, no matter how it's found. Lesson 31 fixes this — not by using a smarter optimizer, but by giving the model more than one projection to work with."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02a906bd",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "1. Increase `lr` in `train` well past a reasonable value (e.g. `lr=20`) on the two-blob dataset. What happens to the loss curve, and how does this relate to the step size overshooting the loss surface's curvature?\n",
    "2. Retrain on the two-blob dataset with several different random seeds. Does the learned direction always end up close to the hand-picked one from Lesson 29, or does it vary a lot? What does that suggest about how many good solutions exist for a well-separated dataset?\n",
    "3. Make the outer ring asymmetric instead of a full circle (e.g. `theta_out = rng.uniform(0, np.pi, 60)`, so the outer class only occupies half the ring). Retrain the single neuron. Does this genuinely broken symmetry change the final accuracy noticeably compared to the ~50% chance-level result on the full ring, and does that support the explanation above (that the *full* ring's symmetric pull on the gradient, not some fundamental inability to learn, is what stalled training near $w \\approx 0$)? (Note: simply translating both classes together by the same offset would *not* break this symmetry, since the bias term $b$ can absorb any shared translation for free.)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.x"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}