Delete lecture_1

lecture_1/README.md

@@ -1,83 +0,0 @@
-# Lecture 1: Introduction and Mathematical Foundations
-
-Chapter 1 of the notes is the toolbox the other six chapters draw on, and it works every piece of
-it on a small example by hand. The six scripts here run those examples. Each one prints the
-chapter's hand-worked value, the value the script computes for it, and the gap between the two, so
-a student reading the chapter can watch its arithmetic come out of the code. Nothing here trains a
-network. The scripts are worked examples, and `tests/test_lecture_1.py` asserts every number they
-print against the chapter that states it.
-
-Run them from the repository root:
-
-```
-cd cis6270
-PYTHONPATH=. python lecture_1/probability.py
-PYTHONPATH=. python -m pytest tests/test_lecture_1.py -q
-```
-
-## The files
-
-| File | Command | What it works through |
-| --- | --- | --- |
-| `probability.py` | `python lecture_1/probability.py` | Sections 1.3 to 1.11. The two-token DNA distribution, the expected adenine count against a Monte Carlo ladder, the worked variance, the three peptides and their covariance matrix, the noising equation, the motif-and-binding table with its marginals and conditionals, Bayes' rule, and the independence check. |
-| `linear_algebra.py` | `python lecture_1/linear_algebra.py` | Sections 1.12, 1.13, 1.15, 1.16, 1.17 and 1.18. The simplex and a mixture, the softmax of `z = (2, 1, 0)`, one step of the two-state chain, the eigenvalues of `A = [[2, 1], [1, 2]]`, powers along an eigendirection, the telescoping `A^n = V Lambda^n V^{-1}`, the stationary distribution and its rate, and `e^{tA}` by both the series and the eigenvalues. |
-| `calculus.py` | `python lecture_1/calculus.py` | Sections 1.22 to 1.30. The Jacobian by automatic differentiation against the chapter's closed form, the determinant as local volume change, change of variables in one and in two dimensions, the log-determinant of a composition, the two inverse branches of the squaring map, the gradient and directional derivative, the divergence of three fields, and the continuity equation checked numerically on `v(x) = x`. |
-| `fields.py` | `python lecture_1/fields.py` | Sections 1.31 and 1.32. The score of the standard Gaussian and of `N(1, 4)`, the normalizing constant differentiated away, the Laplacian of the bowl `x_1^2 + 3 x_2^2` and of the saddle `x_1^2 - x_2^2`, the heat equation at `D = 1/2` on the widening Gaussian `N(0, t)`, and the Fokker-Planck equation with its two reductions, to the continuity equation at `sigma_t = 0` and to the heat equation at `b_t = 0`. The same drift and noise drive an Euler-Maruyama walk, whose sample variance tracks the elapsed time. Chapter 3 starts from both sections. |
-| `information.py` | `python lecture_1/information.py` | Sections 1.33 to 1.35. Information content, entropy, the cross-entropy of `p = (0.8, 0.2)` scored by `q = (0.6, 0.4)`, the split `H(p, q) = H(p) + KL(p \|\| q)`, nonnegativity and asymmetry, and maximum likelihood on ten DNA bases as a KL minimization. |
-| `optimization.py` | `python lecture_1/optimization.py` | Sections 1.36 to 1.38. Gradient descent on `(theta - 3)^2`, an equality constraint through its Lagrangian, normalization as a constraint, the maximum-entropy derivation that produces the exponential family at a GC target of 0.6, the KKT conditions on a bounded quadratic, and the same likelihood problem solved on logits through the softmax. |
-
-## The flags
-
-Every script takes `--steps`, `--seed`, `--out` and `--quiet`, so the flags learned once apply
-across the course. `--steps` drives the iteration count of each script, and every flag a script
-takes changes a line of its report.
-
-| Flag | Where it changes the report |
-| --- | --- |
-| `--steps` | `probability.py`, repeats of the ten-sample Monte Carlo estimate, default 5000. `linear_algebra.py`, iterations of `p_{t+1} = P^T p_t`, default 10. `calculus.py`, Euler steps along the worked flow, default 200. `fields.py`, Euler-Maruyama steps over `[0, 1]`, default 200. `information.py`, repeats of the ten-sample cross-entropy estimate, default 5000. `optimization.py`, gradient-descent steps, default 600. |
-| `--samples` | `probability.py`, the largest Monte Carlo sample size, default 1000. `fields.py`, paths simulated along the SDE, default 20000. `information.py`, the draws behind the sampled cross-entropy and the random distribution pairs, default 20000. `optimization.py`, the competing distributions drawn at the same GC target, default 4000. |
-| `--grid` | `calculus.py` and `fields.py`, the points used in the numerical checks, default 1000. |
-| `--seed` | The Monte Carlo draws, the simulated paths and the random distributions. The same seed gives the same report. |
-| `--out` | A directory for `config.json` and `report.json`, written through `cis6270.runner.save_run`. |
-| `--quiet` | Print the closing report and nothing else. |
-
-## What a default run prints
-
-| Script | Headline numbers |
-| --- | --- |
-| `probability.py` | `E[h(X)] = 1.1000`, Monte Carlo estimate `1.1050` at `N = 1000`, `Var(h(X)) = 0.6900`, spread of 5000 ten-sample averages `0.2616` against the predicted `0.2627`, `rho_XY = 0.9449`, `p(X = 1 \| Y = 1) = 0.8182`, association `0.0960`. |
-| `linear_algebra.py` | `softmax(z)_1 = 0.6652`, `lambda_1 = 3.0000`, `x_10` first entry `118099`, `p*_1 = 0.6000`, second eigenvalue `0.5000`, `e^{tA}_11 = 1.2275` at `t = 0.1`, truncation error `0.0002`. |
-| `calculus.py` | `cos(1) = 0.5403`, `det J_phi = 4.0000` at `x_1 = 2`, `p_Y = 0.1667` inside the image of `A`, `log\|det J\| = 3.1781` for the composition, continuity residual `0.0000`, `log p = -1.0000` along the flow at `t = 1`. |
-| `fields.py` | `s(3) = -3.0000`, `s(5) = -1.0000` under `N(1, 4)`, Laplacian `8.0000` on the bowl and `0.0000` on the saddle, `p_1(0) = 0.3989`, `d p_t(0)/dt = -0.1995`, `D` times the Laplacian `-0.1995`, every residual `0.0000`, sample variance `1.0130` at `t = 1` over 20000 paths. |
-| `information.py` | `H(p) = 0.5004`, `H(p, q) = 0.5919`, `KL(p \|\| q) = 0.0915`, `KL(q \|\| p) = 0.1046`, spread of 5000 ten-sample estimates `0.0502` against the predicted `0.0513`, empirical loss `1.2799` nats, uniform loss `1.3863` nats, `KL(p_hat \|\| uniform) = 0.1064`. |
-| `optimization.py` | `theta = 3.0000` after 600 steps, `lambda = -1.0000` on the equality constraint, `lambda = 0.4055` for a GC target of 0.6, `H(p*) = 1.3662` nats, best drawn competitor `1.3661` nats, `lambda = 4.0000` on the active bound, softmax descent loss `1.2799` nats. |
-
-Every one of these runs in a few seconds on one core. The sampled rows are the ones that move with
-`--seed`. `probability.py --steps 50000` closes the measured spread to `0.2624` and
-`information.py --steps 50000` closes it to `0.0510`, at the cost of a second each.
-
-## What is simplified
-
-- Every script runs in double precision. Chapter 1 quotes five decimal places in places, and the
-  single precision used by the training lectures would not reach them.
-- The Monte Carlo sections redraw their samples, so the sampled values move with `--seed`. The
-  chapter's own ten-sample draw gave 1.2 against the exact 1.1, and that particular draw is not
-  reproduced. The exact value and the `1/sqrt(N)` error scale are.
-- The continuity equation is checked by central differences at interior points of the support of
-  `p_t`, away from the two ends, where a one-sided difference would be needed. The Fokker-Planck
-  residual of `fields.py` uses the same differences, with a second difference for the Laplacian, so
-  it sits near `1e-8` and not at zero.
-- `fields.py` takes the Brownian increment of `dX_t = b_t(X_t) dt + sigma_t dW_t` as a scaled
-  Gaussian draw per step. Chapter 3 defines Brownian motion properly. With a constant drift and a
-  constant noise magnitude the Euler-Maruyama walk has the exact law of the chapter's `N(bt, t)`, so
-  the gap the script prints is Monte Carlo error alone.
-- The branch density `1/(2 sqrt(y))` is integrated numerically over `[0.01, 1]`, because the
-  integrand grows without bound at the origin. Its value over the full `[0, 1]` comes from the
-  antiderivative the chapter gives.
-- The random competitors in the maximum-entropy check and the random pairs in the KL check confirm
-  the chapter's inequalities on a finite sample. They are a check, not the proof, which the chapter
-  gives in full. The competitors cover the two free splits of the GC-constrained family, one inside
-  `{A, T}` and one inside `{C, G}`, drawn independently.
-- Each script has its own three-line `show` helper for the side-by-side printing and its own
-  `WORKED` dictionary of computed values. Neither belongs in `cis6270/`, because no other lecture
-  prints a chapter number beside a computed one.

lecture_1/calculus.py

@@ -1,279 +0,0 @@
-"""Jacobians, change of variables, divergence and the continuity equation.
-
-A map moves points, and the Jacobian is the linear map it looks like near one of them. Its
-determinant is the local volume factor, so the density at an image point is the density at its
-source divided by that factor, the log-determinants add across a composition, and for a
-many-to-one map the density is a sum of one corrected term per inverse branch. In continuous time
-that factor becomes the divergence of a velocity field, and the continuity equation holds it in
-balance with the change in density at every point.
-
-Chapter 1, Sections 1.22 to 1.30 of the notes. Runs on the maps of Sections 1.22.2 and 1.23, the
-squaring map of Section 1.26, and the expanding flow v(x) = x of Section 1.30.2. Every Jacobian
-is taken by automatic differentiation and checked against the chapter's closed form.
-
-    python lecture_1/calculus.py               # every worked example
-    python lecture_1/calculus.py --grid 4000   # finer grids in the numerical checks
-    python lecture_1/calculus.py --steps 2000  # a finer Euler walk along the worked flow
-"""
-
-# %% 01. Imports and configuration
-from __future__ import annotations
-
-import argparse
-import math
-
-import torch
-
-from cis6270.runner import print_report, save_run, set_seed
-
-torch.set_default_dtype(torch.float64)  # Chapter 1 is arithmetic, and five printed digits need it.
-
-
-A_SCALE = torch.tensor([[2.0, 0.0], [0.0, 3.0]])  # The map of Section 1.23, which scales area by 6.
-
-
-WORKED: dict[str, float] = {}  # Every computed value, by label, for tests/test_lecture_1.py.
-
-
-def show(label: str, by_hand: float, computed: float, quiet: bool = False) -> None:
-    """Print the chapter's hand-worked value beside the value this file computes, and keep it."""
-    WORKED[label] = computed  # The test suite asserts these against the chapter's own numbers.
-    if not quiet:  # Under --quiet only the closing report is printed.
-        computed = computed + 0.0  # Adding zero prints a negative zero as a plain zero.
-        print(f"  {label:<36s} chapter {by_hand:>11.5f}   computed {computed:>11.5f}"
-              f"   gap {abs(by_hand - computed):.1e}")
-
-
-def phi_mixed(x: torch.Tensor) -> torch.Tensor:
-    """The map (x_1^2 + x_2, sin x_1 + 3 x_2) of Section 1.22.2, nonlinear in its first component."""
-    return torch.stack([x[0] ** 2 + x[1], torch.sin(x[0]) + 3.0 * x[1]])
-
-
-def phi_square(x: torch.Tensor) -> torch.Tensor:
-    """The map (x_1^2, x_2) of Section 1.23, one-to-one once the domain is restricted to x_1 > 0."""
-    return torch.stack([x[0] ** 2, x[1]])
-
-
-def psi_scale(y: torch.Tensor) -> torch.Tensor:
-    """The linear map y -> A y with A = diag(2, 3), composed after phi_square in Section 1.25."""
-    return A_SCALE @ y
-
-
-def f_bowl(x: torch.Tensor) -> torch.Tensor:
-    """The scalar field f(x_1, x_2) = x_1^2 + 2 x_2^2 of Section 1.27.1."""
-    return x[0] ** 2 + 2.0 * x[1] ** 2
-
-
-def divergence(field, x: torch.Tensor) -> float:
-    """Section 1.29: the divergence is the sum of d v_i / d x_i, which is the trace of J_v."""
-    return float(torch.autograd.functional.jacobian(field, x).trace())
-
-
-# %% 02. The Jacobian matrix and the local linear approximation
-def jacobian(quiet: bool = False) -> dict:
-    """Sections 1.22.2 and 1.22.3: J_phi at (1, 0), and the displacement it predicts."""
-    point = torch.tensor([1.0, 0.0])  # The point at which the chapter evaluates the Jacobian.
-    computed = torch.autograd.functional.jacobian(phi_mixed, point)  # [2, 2], by backpropagation.
-    closed_form = torch.tensor([[2.0 * point[0], 1.0],  # The chapter's entries, 2 x_1 and 1.
-                                [torch.cos(point[0]), 3.0]])  # cos x_1 and 3; [2, 2].
-    delta = torch.tensor([0.01, -0.02])  # The small input displacement of Section 1.22.3.
-    predicted = computed @ delta  # J_phi(x) Delta x, the first-order change in the output; [2].
-    exact = phi_mixed(point + delta) - phi_mixed(point)  # The change the full map produces; [2].
-    show("J_phi(1, 0)_11", 2.0, float(computed[0, 0]), quiet)
-    show("J_phi(1, 0)_21 = cos(1)", 0.5403, float(computed[1, 0]), quiet)
-    show("largest |autograd - closed form|", 0.0, float((computed - closed_form).abs().max()), quiet)
-    show("predicted change, first output", 0.0, float(predicted[0]), quiet)
-    show("predicted change, second output", -0.0546, float(predicted[1]), quiet)
-    show("exact change, first output", 0.0001, float(exact[0]), quiet)
-    show("exact change, second output", -0.05464, float(exact[1]), quiet)
-    return {"jacobian": computed, "predicted": predicted, "exact": exact}
-
-
-# %% 03. The Jacobian determinant as local volume change
-def determinant(quiet: bool = False) -> dict:
-    """Section 1.23: the absolute determinant is the factor by which a small region changes area."""
-    corners = torch.tensor([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])  # The unit square.
-    mapped = corners @ A_SCALE.T  # The image rectangle [0, 2] by [0, 3]; [4, 2].
-    area = float((mapped[1, 0] - mapped[0, 0]) * (mapped[3, 1] - mapped[0, 1]))  # Width times height.
-    collapsing = torch.tensor([[1.0, 0.0], [0.0, 0.0]])  # A map that flattens the plane onto a line.
-    stretch = torch.autograd.functional.jacobian(phi_square, torch.tensor([2.0, 1.0]))  # At x_1 = 2.
-    small = torch.tensor([0.25, 1.0])  # A point where the same map contracts, since 2 x_1 is below 1.
-    squeeze = torch.autograd.functional.jacobian(phi_square, small)  # J_phi at x_1 = 0.25; [2, 2].
-    show("det A for diag(2, 3)", 6.0, float(torch.linalg.det(A_SCALE)), quiet)
-    show("area of the image rectangle", 6.0, area, quiet)
-    show("det of the collapsing map", 0.0, float(torch.linalg.det(collapsing)), quiet)
-    show("det J_phi at x_1 = 2", 4.0, float(torch.linalg.det(stretch)), quiet)
-    show("det J_phi at x_1 = 0.25", 0.5, float(torch.linalg.det(squeeze)), quiet)
-    return {"area": area, "stretch": float(torch.linalg.det(stretch))}
-
-
-# %% 04. Change of variables in one and in two dimensions
-def change_of_variables(grid: int, quiet: bool = False) -> dict:
-    """Sections 1.24.1 and 1.24.2: p_Y(phi(x)) = p_X(x) / |det J_phi(x)| on the two worked maps."""
-    inverse_slope = 0.5  # d phi^{-1}/dy for phi^{-1}(y) = (y - 1)/2, the reciprocal of phi' = 2.
-    p_y = 1.0 * abs(inverse_slope)  # p_Y(y) = p_X(phi^{-1}(y)) |d phi^{-1}/dy|, so p_Y = 1/2.
-    preimage = ((1.4 - 1.0) / 2.0, (2.0 - 1.0) / 2.0)  # The inverse image of B = [1.4, 2.0].
-    mass = (preimage[1] - preimage[0]) * 1.0  # P(X in [0.2, 0.5]) under the uniform density.
-    line = 1.0 + (torch.arange(grid) + 0.5) * (2.0 / grid)  # Midpoints of [1, 3]; [grid].
-    normalization = float(torch.full_like(line, p_y).sum() * (2.0 / grid))  # Integral of p_Y there.
-    flat = 1.0 / float(torch.linalg.det(A_SCALE))  # p_Y = 1/6 inside the image of the unit square.
-    step = 1.0 / grid  # Midpoint spacing on the unit interval of the first input coordinate.
-    x_one = (torch.arange(grid) + 0.5) * step  # Midpoints of (0, 1) in that coordinate; [grid].
-    volume_factor = 2.0 * x_one  # |det J_phi| = 2 x_1 for phi = (x_1^2, x_2), at each point; [grid].
-    output_density = 1.0 / (2.0 * (x_one ** 2).sqrt())  # p_Y(y) = 1/(2 sqrt(y_1)) at y_1 = x_1^2.
-    residual = float((output_density * volume_factor - 1.0).abs().max())  # Should recover p_X = 1.
-    transported = float((output_density * volume_factor).sum() * step)  # Total mass by substitution.
-    show("preimage of [1.4, 2.0], lower end", 0.2, preimage[0], quiet)
-    show("preimage of [1.4, 2.0], upper end", 0.5, preimage[1], quiet)
-    show("P(X in the preimage)", 0.3, mass, quiet)
-    show("p_Y on [1, 3] after Y = 2X + 1", 0.5, p_y, quiet)
-    show("integral of p_Y over [1, 3]", 1.0, normalization, quiet)
-    show("p_Y inside the image of A", 1.0 / 6.0, flat, quiet)
-    show("density times area for A", 1.0, flat * 6.0, quiet)
-    show("p_Y(y_1 = 1) for the x_1^2 map", 0.5, 1.0 / (2.0 * math.sqrt(1.0)), quiet)
-    show("largest |p_Y |det J| - p_X|", 0.0, residual, quiet)
-    show("total mass after the x_1^2 map", 1.0, transported, quiet)
-    return {"p_y": p_y, "flat": flat, "residual": residual, "transported": transported}
-
-
-# %% 05. Composed maps and the log-density form
-def composition(quiet: bool = False) -> dict:
-    """Section 1.25: log|det J| adds across a composition, so two factors become one sum."""
-    point = torch.tensor([2.0, 1.0])  # The chapter evaluates the composition at x_1 = 2.
-    inner = torch.autograd.functional.jacobian(phi_square, point)  # J_phi(x), det 2 x_1 = 4; [2, 2].
-    outer = torch.autograd.functional.jacobian(psi_scale, phi_square(point))  # J_psi at phi(x); [2, 2].
-    chained = outer @ inner  # The chain rule of Proposition 1.1, giving diag(8, 3); [2, 2].
-    composed = torch.autograd.functional.jacobian(lambda x: psi_scale(phi_square(x)), point)
-    log_inner = math.log(abs(float(torch.linalg.det(inner))))  # log 4 = 1.3863.
-    log_outer = math.log(abs(float(torch.linalg.det(outer))))  # log 6 = 1.7918.
-    log_density = math.log(1.0) - (log_inner + log_outer)  # log p_Y = log p_X - log|det J|.
-    show("det J_psi", 6.0, float(torch.linalg.det(outer)), quiet)
-    show("chained Jacobian, entry 11", 8.0, float(chained[0, 0]), quiet)
-    show("det of the composition", 24.0, float(torch.linalg.det(composed)), quiet)
-    show("log 4", 1.3863, log_inner, quiet)
-    show("log 6", 1.7918, log_outer, quiet)
-    show("log 4 + log 6 = log 24", 3.1781, log_inner + log_outer, quiet)
-    show("log p_Y from a unit base density", -3.1781, log_density, quiet)
-    return {"log_sum": log_inner + log_outer, "log_density": log_density}
-
-
-# %% 06. Local invertibility and several inverse branches
-def branches(grid: int, quiet: bool = False) -> dict:
-    """Section 1.26: for Y = X^2 on a uniform X, the two preimages each contribute to p_Y(y)."""
-    def branch_density(y: float) -> float:
-        roots = (math.sqrt(y), -math.sqrt(y))  # The two solutions of x^2 = y inside the support.
-        return sum(0.5 / abs(2.0 * root) for root in roots)  # Sum of p_X(x) / |phi'(x)| over roots.
-
-    lower = 0.01  # The integral is cut at 0.01 because the integrand grows without bound at zero.
-    step = (1.0 - lower) / grid  # Midpoint spacing on [0.01, 1].
-    points = lower + (torch.arange(grid) + 0.5) * step  # Midpoints of the integration interval.
-    partial = float((1.0 / (2.0 * points.sqrt())).sum() * step)  # Midpoint integral of 1/(2 sqrt y).
-    show("p_Y(0.25), two branches at ±0.5", 1.0, branch_density(0.25), quiet)
-    show("p_Y(0.04), two branches at ±0.2", 2.5, branch_density(0.04), quiet)
-    show("integral of p_Y over [0.01, 1]", 1.0 - math.sqrt(lower), partial, quiet)
-    return {"at_quarter": branch_density(0.25), "partial": partial}
-
-
-# %% 07. Scalar fields, the gradient and the directional derivative
-def gradient(quiet: bool = False) -> dict:
-    """Sections 1.27.1 and 1.27.2: grad f at (1, 2), its steepest direction and a directional rate."""
-    point = torch.tensor([1.0, 2.0])  # The point at which the chapter evaluates the gradient.
-    grad = torch.autograd.functional.jacobian(f_bowl, point)  # (2 x_1, 4 x_2) = (2, 8); [2].
-    steepest = grad / grad.norm()  # The unit vector of fastest increase; [2].
-    along_x1 = float(grad @ torch.tensor([1.0, 0.0]))  # D_u f for u the first basis direction.
-    along_x2 = float(grad @ torch.tensor([0.0, 1.0]))  # D_u f for u the second basis direction.
-    show("grad f(1, 2), first entry", 2.0, float(grad[0]), quiet)
-    show("grad f(1, 2), second entry", 8.0, float(grad[1]), quiet)
-    show("|grad f(1, 2)| = sqrt(68)", math.sqrt(68.0), float(grad.norm()), quiet)
-    show("steepest direction, first entry", 0.243, float(steepest[0]), quiet)
-    show("steepest direction, second entry", 0.970, float(steepest[1]), quiet)
-    show("D_u f along (1, 0)", 2.0, along_x1, quiet)
-    show("D_u f along (0, 1)", 8.0, along_x2, quiet)
-    return {"gradient": grad, "along_x1": along_x1}
-
-
-# %% 08. Divergence
-def divergences(quiet: bool = False) -> dict:
-    """Section 1.29: the divergence of the chapter's three fields is -3, then 1, then 0."""
-    point = torch.tensor([1.0, 1.0])  # The divergences below are the same at every point.
-    contract = divergence(lambda x: torch.stack([-x[0], -2.0 * x[1]]), point)  # -1 + (-2) = -3.
-    expand = divergence(lambda x: torch.stack([2.0 * x[0], -x[1]]), point)  # 2 - 1 = 1.
-    rotate = divergence(lambda x: torch.stack([-x[1], x[0]]), point)  # 0 + 0 = 0.
-    show("div of (-x_1, -2 x_2)", -3.0, contract, quiet)
-    show("div of (2 x_1, -x_2)", 1.0, expand, quiet)
-    show("div of (-x_2, x_1)", 0.0, rotate, quiet)
-    show("area factor e^{div t} at t = 0.5", 1.649, math.exp(expand * 0.5), quiet)
-    return {"contract": contract, "expand": expand, "rotate": rotate}
-
-
-# %% 09. The continuity equation on the worked flow
-def continuity(grid: int, steps: int, quiet: bool = False) -> dict:
-    """Section 1.30.2: v(x) = x stretches a uniform source, with p_t(x) = e^{-t} on [0, e^t].
-
-    Central differences give the continuity residual at interior points, and one Euler walk
-    accumulates the log-density rate -div v of Section 1.30.1 along a trajectory.
-    """
-    def flux(time: float, position):
-        return math.exp(-time) * position  # j_t(x) = p_t(x) v_t(x), with p_t = e^{-t} and v(x) = x.
-
-    t, h = 0.5, 1e-4  # The time of the residual check, and the central-difference step.
-    x = torch.linspace(0.1, 1.4, grid)  # Interior points of the support [0, e^{0.5}]; [grid].
-    time_rate = (math.exp(-(t + h)) - math.exp(-(t - h))) / (2.0 * h)  # d p_t / d t at fixed x.
-    flux_rate = (flux(t, x + h) - flux(t, x - h)) / (2.0 * h)  # d (p_t v) / d x; [grid].
-    residual = float((time_rate + flux_rate).abs().max())  # The continuity equation, point by point.
-    x_t, log_p, dt = 0.5, 0.0, 1.0 / max(steps, 1)  # One trajectory started inside the source.
-    for _ in range(max(steps, 1)):
-        log_p -= 1.0 * dt  # d/dt log p_t(x_t) = -div v = -1 for the field v(x) = x.
-        x_t = x_t + x_t * dt  # Euler step on dx_t/dt = v(x_t) = x_t.
-    show("largest continuity residual", 0.0, residual, quiet)
-    show("log p along the trajectory at t = 1", -1.0, log_p, quiet)
-    show("p_t(x_t) at t = 1", 0.3679, math.exp(log_p), quiet)
-    show("x_t at t = 1, exact 0.5e", 0.5 * math.e, x_t, quiet)
-    show("density at t = log 2", 0.5, math.exp(-math.log(2.0)), quiet)
-    show("mass at t = log 2", 1.0, 2.0 * math.exp(-math.log(2.0)), quiet)
-    show("mass at t = 1", 1.000, math.e * math.exp(-1.0), quiet)
-    return {"residual": residual, "log_p": log_p, "x_t": x_t}
-
-
-# %% 10. The report
-def main() -> None:
-    parser = argparse.ArgumentParser(description=__doc__.split("\n")[0])
-    parser.add_argument("--grid", type=int, default=1000, help="points in the numerical checks")
-    parser.add_argument("--steps", type=int, default=200, help="Euler steps along the worked flow")
-    parser.add_argument("--seed", type=int, default=0, help="random seed, for a repeatable run")
-    parser.add_argument("--out", type=str, default="", help="directory for the run's output")
-    parser.add_argument("--quiet", action="store_true", help="print only the final report")
-    args = parser.parse_args()
-
-    set_seed(args.seed)  # Seed first, so the same seed gives the same run.
-    if not args.quiet:
-        print("Chapter 1, Sections 1.22 to 1.30, worked example by worked example\n")
-    local = jacobian(args.quiet)
-    volume = determinant(args.quiet)
-    moved = change_of_variables(args.grid, args.quiet)
-    stacked = composition(args.quiet)
-    folded = branches(args.grid, args.quiet)
-    slope = gradient(args.quiet)
-    fields = divergences(args.quiet)
-    flow = continuity(args.grid, args.steps, args.quiet)
-
-    report = {
-        "J_phi(1, 0)_21 = cos(1)": float(local["jacobian"][1, 0]),
-        "det J_phi at x_1 = 2": volume["stretch"],
-        "p_Y on [1, 3] after Y = 2X + 1": moved["p_y"],
-        "p_Y inside the image of A": moved["flat"],
-        "mass after the x_1^2 map": moved["transported"],
-        "log|det J| of the composition": stacked["log_sum"],
-        "p_Y(0.25) summed over branches": folded["at_quarter"],
-        "D_u f along (1, 0)": slope["along_x1"],
-        "div of (2 x_1, -x_2)": fields["expand"],
-        "largest continuity residual": flow["residual"],
-        "log p along the flow at t = 1": flow["log_p"],
-        "x_t at t = 1 after Euler steps": flow["x_t"],
-    }
-    print_report("Calculus, Sections 1.22 to 1.30", report)
-    save_run(args.out, config=vars(args), report=report)
-
-
-if __name__ == "__main__":  # Execute when launched directly.
-    main()  # Work every example of Sections 1.22 to 1.30 and print the report.

lecture_1/fields.py

@@ -1,234 +0,0 @@
-"""Score functions, the Laplacian and the Fokker-Planck equation.
-
-The score of a density is the gradient of its logarithm, a vector field of the same shape as a
-velocity field, and the normalizing constant differentiates away. The Laplacian is the divergence
-of the gradient, and it measures the curvature that drives diffusion, so the heat equation sets the
-rate of change of a density proportional to it. One statement about the density holds both
-mechanisms at once, the Fokker-Planck equation, and it returns the continuity equation when the
-noise magnitude is zero and the heat equation when the drift is.
-
-Chapter 1, Sections 1.31 and 1.32 of the notes. Runs on the Gaussians of Section 1.31, the bowl and
-the saddle of Section 1.32, the widening Gaussian p_t = N(0, t) of Section 1.32.1 and the expanding
-flow v(x) = x of Section 1.30.2. Prints the chapter's number beside the computed one.
-
-    python lecture_1/fields.py               # every worked example
-    python lecture_1/fields.py --steps 2000  # a finer Euler-Maruyama walk
-    python lecture_1/fields.py --quiet       # the closing report only
-"""
-
-# %% 01. Imports and configuration
-from __future__ import annotations
-
-import argparse
-import math
-
-import torch
-
-from cis6270.runner import print_report, save_run, set_seed
-
-torch.set_default_dtype(torch.float64)  # Chapter 1 is arithmetic, and five printed digits need it.
-
-
-SCORE_POINTS = torch.tensor([3.0, -2.0, 0.0])  # The three states scored in Section 1.31.
-WIDE_MEAN, WIDE_VARIANCE = 1.0, 4.0  # The general Gaussian N(mu, sigma^2) of Section 1.31.
-DRIFT = 0.5  # The constant drift b_t(x) = 0.5 of the combined process in Section 1.32.2.
-NOISE = 1.0  # The noise magnitude sigma_t, so that D = sigma_t^2 / 2 = 0.5 as in Section 1.32.1.
-
-
-WORKED: dict[str, float] = {}  # Every computed value, by label, for tests/test_lecture_1.py.
-
-
-def show(label: str, by_hand: float, computed: float, quiet: bool = False) -> None:
-    """Print the chapter's hand-worked value beside the value this file computes, and keep it."""
-    WORKED[label] = computed  # The test suite asserts these against the chapter's own numbers.
-    if not quiet:  # Under --quiet only the closing report is printed.
-        computed = computed + 0.0  # Adding zero prints a negative zero as a plain zero.
-        print(f"  {label:<36s} chapter {by_hand:>11.5f}   computed {computed:>11.5f}"
-              f"   gap {abs(by_hand - computed):.1e}")
-
-
-def log_gaussian(x: torch.Tensor, mu: float = 0.0, variance: float = 1.0) -> torch.Tensor:
-    """log N(x; mu, sigma^2) = -(x - mu)^2 / (2 sigma^2) - (1/2) log(2 pi sigma^2)."""
-    normalizer = 0.5 * math.log(2.0 * math.pi * variance)  # The term with no x in it.
-    return -((x - mu) ** 2) / (2.0 * variance) - normalizer
-
-
-def gaussian(t: float, x: torch.Tensor) -> torch.Tensor:
-    """p_t = N(0, t) in one dimension, the widening Gaussian of Section 1.32.1."""
-    return (-(x ** 2) / (2.0 * t)).exp() / math.sqrt(2.0 * math.pi * t)
-
-
-def score(log_density, x: torch.Tensor) -> torch.Tensor:
-    """Section 1.31: s(x) = grad_x log p(x), taken by automatic differentiation."""
-    point = x.detach().clone().requires_grad_(True)  # A leaf to differentiate against; [N].
-    total = log_density(point).sum()  # The points do not interact, so one pass gives every score.
-    gradient, = torch.autograd.grad(total, point)  # d/dx log p at each point; [N].
-    return gradient
-
-
-def laplacian(field, x: torch.Tensor) -> float:
-    """Section 1.32: the Laplacian is the divergence of the gradient, the trace of the Hessian."""
-    return float(torch.autograd.functional.hessian(field, x).trace())
-
-
-def fokker_planck_residual(p, b, sigma: float, x: torch.Tensor, t: float) -> float:
-    """Section 1.32.2: d p_t / d t + div(b_t p_t) - (1/2) sigma_t^2 Laplacian p_t, by differences.
-
-    Setting sigma to zero leaves the continuity equation of Section 1.30, and setting b to zero
-    leaves the heat equation of Section 1.32.1, so one residual covers all three statements.
-    """
-    h = 1e-4  # The central-difference step, small enough for five digits in double precision.
-    time_rate = (p(t + h, x) - p(t - h, x)) / (2.0 * h)  # d p_t(x) / d t at fixed x; [grid].
-    flux = (b(x + h) * p(t, x + h) - b(x - h) * p(t, x - h)) / (2.0 * h)  # div(b_t p_t); [grid].
-    curvature = (p(t, x + h) - 2.0 * p(t, x) + p(t, x - h)) / h ** 2  # Laplacian p_t; [grid].
-    return float((time_rate + flux - 0.5 * sigma ** 2 * curvature).abs().max())
-
-
-# %% 02. The score of a Gaussian
-def scores(quiet: bool = False) -> dict:
-    """Section 1.31: the score of N(0, 1) is -x, and of N(mu, sigma^2) it is -(x - mu)/sigma^2."""
-    standard = score(log_gaussian, SCORE_POINTS)  # s(x) at x = 3, -2 and 0; [3].
-    closed_form = -SCORE_POINTS  # The chapter's result, with the log(2 pi) term differentiated away.
-    unnormalized = score(lambda z: -(z ** 2) / 2.0, SCORE_POINTS)  # The same density without its Z.
-    wide = score(lambda z: log_gaussian(z, WIDE_MEAN, WIDE_VARIANCE), torch.tensor([5.0]))  # [1].
-    show("s(x) at x = 3", -3.0, float(standard[0]), quiet)
-    show("s(x) at x = -2", 2.0, float(standard[1]), quiet)
-    show("s(x) at the mode x = 0", 0.0, float(standard[2]), quiet)
-    show("largest |s(x) + x|", 0.0, float((standard - closed_form).abs().max()), quiet)
-    show("largest change from dropping Z", 0.0, float((unnormalized - standard).abs().max()), quiet)
-    show("s(5) under N(1, 4)", -1.0, float(wide[0]), quiet)
-    return {"standard": standard, "wide": float(wide[0])}
-
-
-# %% 03. The Laplacian of a scalar field
-def curvatures(quiet: bool = False) -> dict:
-    """Section 1.32: the bowl x_1^2 + 3 x_2^2 has Laplacian 8, and the saddle x_1^2 - x_2^2 has 0."""
-    def bowl(x: torch.Tensor) -> torch.Tensor:
-        return x[0] ** 2 + 3.0 * x[1] ** 2  # The field whose two curvatures are 2 and 6.
-
-    def saddle(x: torch.Tensor) -> torch.Tensor:
-        return x[0] ** 2 - x[1] ** 2  # The field whose two curvatures cancel.
-
-    point = torch.tensor([1.0, 2.0])  # The Laplacians below are the same at every point.
-    second = torch.autograd.functional.hessian(bowl, point)  # The two pure second derivatives; [2, 2].
-    h = 1e-3  # The spacing of the four neighbors used for the average reading of the Laplacian.
-    steps = torch.tensor([[h, 0.0], [-h, 0.0], [0.0, h], [0.0, -h]])  # The four offsets; [4, 2].
-    neighbors = sum(bowl(point + step) for step in steps) - 4.0 * bowl(point)  # Above the center.
-    show("d^2 f / d x_1^2", 2.0, float(second[0, 0]), quiet)
-    show("d^2 f / d x_2^2", 6.0, float(second[1, 1]), quiet)
-    show("Laplacian of x_1^2 + 3 x_2^2", 8.0, laplacian(bowl, point), quiet)
-    show("Laplacian from the four neighbors", 8.0, float(neighbors) / h ** 2, quiet)
-    show("Laplacian of the saddle", 0.0, laplacian(saddle, point), quiet)
-    return {"bowl": laplacian(bowl, point), "saddle": laplacian(saddle, point)}
-
-
-# %% 04. The heat equation
-def heat(grid: int, quiet: bool = False) -> dict:
-    """Section 1.32.1: at D = 1/2 the Gaussian N(0, t) solves d p_t / d t = D Laplacian p_t."""
-    diffusion = 0.5 * NOISE ** 2  # D = sigma_t^2 / 2, the coefficient the chapter takes as 1/2.
-    center = 1.0 / math.sqrt(2.0 * math.pi)  # p_1(0), the height at the middle at t = 1.
-    time_rate = -0.5 / math.sqrt(2.0 * math.pi)  # d p_t(0)/dt = -(1/2)(2 pi)^{-1/2} t^{-3/2} at t = 1.
-    curvature = laplacian(lambda z: gaussian(1.0, z), torch.tensor([0.0]))  # Laplacian p_1 at 0.
-    x = torch.linspace(-3.0, 3.0, grid)  # Points covering the bulk of the mass at t = 1; [grid].
-    residual = fokker_planck_residual(
-        gaussian, torch.zeros_like, NOISE, x, 1.0)  # The heat equation, with no drift.
-    show("p_1(0) = 1 / sqrt(2 pi)", 0.39894, center, quiet)
-    show("d p_t(0) / d t at t = 1", -0.19947, time_rate, quiet)
-    show("Laplacian of p_1 at 0", -0.39894, curvature, quiet)
-    show("-p_1(0) / sigma^2 at sigma^2 = 1", -0.39894, -center / 1.0, quiet)
-    show("D times the Laplacian of p_1 at 0", -0.19947, diffusion * curvature, quiet)
-    show("largest heat-equation residual", 0.0, residual, quiet)
-    return {"center": center, "time_rate": time_rate, "curvature": curvature,
-            "flow": diffusion * curvature, "residual": residual}
-
-
-# %% 05. Drift and diffusion together
-def fokker_planck(grid: int, quiet: bool = False) -> dict:
-    """Section 1.32.2: a drift that translates and a noise that widens, in one equation.
-
-    With a constant drift b and sigma_t = 1 the density N(bt, t) solves the equation, since the
-    drift moves the center and the noise grows the variance at the rate of the heat equation.
-    """
-    def translating(t: float, x: torch.Tensor) -> torch.Tensor:
-        return gaussian(t, x - DRIFT * t)  # N(bt, t), the center carried along by the drift.
-
-    def expanding(t: float, x: torch.Tensor) -> torch.Tensor:
-        return torch.full_like(x, math.exp(-t))  # p_t = e^{-t} on [0, e^t], Section 1.30.2.
-
-    x = torch.linspace(-3.0, 3.0, grid)  # Points across the translating Gaussian; [grid].
-    both = fokker_planck_residual(translating, lambda z: torch.full_like(z, DRIFT), NOISE, x, 1.0)
-    inside = torch.linspace(0.1, 1.4, grid)  # Interior points of the support at t = 0.5; [grid].
-    transport = fokker_planck_residual(expanding, lambda z: z, 0.0, inside, 0.5)  # sigma_t = 0.
-    center = DRIFT * 1.0  # The mean of p_t at t = 1, which the drift alone sets.
-    show("center of p_t at t = 1", 0.5, center, quiet)
-    show("largest residual, drift and noise", 0.0, both, quiet)
-    show("largest residual at sigma_t = 0", 0.0, transport, quiet)
-    return {"both": both, "transport": transport, "center": center}
-
-
-# %% 06. The process behind the density
-def euler_maruyama(samples: int, steps: int, quiet: bool = False) -> dict:
-    """Section 1.32.2: dX_t = b_t(X_t) dt + sigma_t dW_t, simulated, against the density above.
-
-    The variance of p_t = N(bt, t) equals the elapsed time, so the sample mean and the sample
-    variance of the walk track bt and t. Both gaps shrink as --samples grows.
-    """
-    dt = 1.0 / max(steps, 1)  # The time increment of one Euler-Maruyama step.
-    x_t = torch.zeros(max(samples, 1))  # Every path starts at the origin, so p_0 is a point mass.
-    mean_gap, variance_gap = 0.0, 0.0  # The largest departure from bt and from t along the walk.
-    for step in range(max(steps, 1)):
-        noise = torch.randn_like(x_t)  # The Brownian increment dW_t, scaled by sqrt(dt); [samples].
-        x_t = x_t + DRIFT * dt + NOISE * math.sqrt(dt) * noise  # One step of the chapter's SDE.
-        elapsed = (step + 1) * dt  # The time reached after this step.
-        mean_gap = max(mean_gap, abs(float(x_t.mean()) - DRIFT * elapsed))
-        variance_gap = max(variance_gap, abs(float(x_t.var(unbiased=False)) - elapsed))
-    show("sample mean at t = 1, exact 0.5", 0.5, float(x_t.mean()), quiet)
-    show("sample variance at t = 1, exact 1", 1.0, float(x_t.var(unbiased=False)), quiet)
-    show("largest |sample mean - b t|", 0.0, mean_gap, quiet)
-    show("largest |sample variance - t|", 0.0, variance_gap, quiet)
-    return {"mean": float(x_t.mean()), "variance": float(x_t.var(unbiased=False)),
-            "mean_gap": mean_gap, "variance_gap": variance_gap, "steps": max(steps, 1)}
-
-
-# %% 07. The report
-def main() -> None:
-    parser = argparse.ArgumentParser(description=__doc__.split("\n")[0])
-    parser.add_argument("--grid", type=int, default=1000, help="points in the numerical checks")
-    parser.add_argument("--steps", type=int, default=200, help="Euler-Maruyama steps over [0, 1]")
-    parser.add_argument("--samples", type=int, default=20000, help="paths simulated along the SDE")
-    parser.add_argument("--seed", type=int, default=0, help="random seed for the simulated paths")
-    parser.add_argument("--out", type=str, default="", help="directory for the run's output")
-    parser.add_argument("--quiet", action="store_true", help="print only the final report")
-    args = parser.parse_args()
-
-    set_seed(args.seed)  # Seed first, so the same seed gives the same walk.
-    if not args.quiet:
-        print("Chapter 1, Sections 1.31 and 1.32, worked example by worked example\n")
-    uphill = scores(args.quiet)
-    bowl = curvatures(args.quiet)
-    spreading = heat(args.grid, args.quiet)
-    combined = fokker_planck(args.grid, args.quiet)
-    walk = euler_maruyama(args.samples, args.steps, args.quiet)
-
-    report = {
-        "s(x) at x = 3": float(uphill["standard"][0]),
-        "s(5) under N(1, 4)": uphill["wide"],
-        "Laplacian of x_1^2 + 3 x_2^2": bowl["bowl"],
-        "Laplacian of the saddle": bowl["saddle"],
-        "p_1(0)": spreading["center"],
-        "d p_t(0) / d t at t = 1": spreading["time_rate"],
-        "Laplacian of p_1 at 0": spreading["curvature"],
-        "D times that Laplacian": spreading["flow"],
-        "heat-equation residual": spreading["residual"],
-        "Fokker-Planck residual": combined["both"],
-        "residual at sigma_t = 0": combined["transport"],
-        "Euler-Maruyama steps": walk["steps"],
-        "sample variance at t = 1": walk["variance"],
-        "largest |sample variance - t|": walk["variance_gap"],
-    }
-    print_report("Fields, Sections 1.31 and 1.32", report)
-    save_run(args.out, config=vars(args), report=report)
-
-
-if __name__ == "__main__":  # Execute when launched directly.
-    main()  # Work every example of Sections 1.31 and 1.32 and print the report.

lecture_1/information.py

@@ -1,219 +0,0 @@
-"""Entropy, cross-entropy, KL divergence and maximum likelihood, in nats.
-
-Information content measures one observation, entropy averages it under the distribution that
-produced it, cross-entropy averages it under a second distribution supplying the probabilities,
-and the KL divergence is the difference between those two averages. Minimizing a negative
-log-likelihood over a model is therefore minimizing a KL divergence to the data distribution.
-
-Chapter 1, Sections 1.33, 1.34 and 1.35 of the notes. Runs on the pair p = (0.8, 0.2) and
-q = (0.6, 0.4) of Section 1.33.3 and on the ten DNA bases of Section 1.35. Prints the chapter's
-number beside the computed one. Logarithms are natural, so every quantity is in nats.
-
-    python lecture_1/information.py                   # every worked example
-    python lecture_1/information.py --samples 200000  # a tighter Monte Carlo cross-entropy
-    python lecture_1/information.py --quiet           # the closing report only
-"""
-
-# %% 01. Imports and configuration
-from __future__ import annotations
-
-import argparse
-import math
-
-import torch
-
-from cis6270.runner import print_report, save_run, set_seed
-
-torch.set_default_dtype(torch.float64)  # Chapter 1 is arithmetic, and five printed digits need it.
-
-
-P_REF = torch.tensor([0.8, 0.2])  # The reference distribution p of Sections 1.33.3 and 1.34.
-Q_MODEL = torch.tensor([0.6, 0.4])  # The comparison distribution q scored against those samples.
-BASE_COUNTS = torch.tensor([4.0, 3.0, 2.0, 1.0])  # n(A), n(C), n(G), n(T) over ten DNA bases.
-
-
-WORKED: dict[str, float] = {}  # Every computed value, by label, for tests/test_lecture_1.py.
-
-
-def show(label: str, by_hand: float, computed: float, quiet: bool = False) -> None:
-    """Print the chapter's hand-worked value beside the value this file computes, and keep it."""
-    WORKED[label] = computed  # The test suite asserts these against the chapter's own numbers.
-    if not quiet:  # Under --quiet only the closing report is printed.
-        computed = computed + 0.0  # Adding zero prints a negative zero as a plain zero.
-        print(f"  {label:<36s} chapter {by_hand:>11.5f}   computed {computed:>11.5f}"
-              f"   gap {abs(by_hand - computed):.1e}")
-
-
-def entropy(p: torch.Tensor) -> float:
-    """H(p) = -sum_x p(x) log p(x), the expected information in one sample drawn from p."""
-    return float(-(p * p.clamp_min(1e-300).log()).sum())  # Clamp keeps 0 log 0 at zero.
-
-
-def cross_entropy(p: torch.Tensor, q: torch.Tensor) -> float:
-    """H(p, q) = -sum_x p(x) log q(x), the cost of scoring samples from p with the values of q."""
-    return float(-(p * q.log()).sum())  # The weights come from p and the probabilities from q.
-
-
-def divergence(p: torch.Tensor, q: torch.Tensor) -> float:
-    """KL(p || q) = sum_x p(x) log(p(x) / q(x)), the expected log ratio taken under p."""
-    return float((p * (p / q).log()).sum())  # Zero exactly when q matches p at every state.
-
-
-# %% 02. Information content
-def information(quiet: bool = False) -> dict:
-    """Section 1.33.1: I(x) = -log p(x), so a rarer outcome carries more information."""
-    common = -math.log(0.5)  # An outcome of probability one half.
-    rare = -math.log(0.01)  # An outcome fifty times less likely.
-    certain = -math.log(1.0)  # A certain event, which has zero information.
-    show("I(x) at p(x) = 0.5", 0.693, common, quiet)
-    show("I(x) at p(x) = 0.01", 4.605, rare, quiet)
-    show("I(x) at p(x) = 1", 0.0, certain, quiet)
-    show("ratio of the two surprises", 6.6, rare / common, quiet)
-    return {"common": common, "rare": rare}
-
-
-# %% 03. Entropy
-def entropies(quiet: bool = False) -> dict:
-    """Section 1.33.2: entropy is high when the mass is spread and zero when it is all on one state."""
-    skewed = entropy(torch.tensor([0.9, 0.1]))  # 0.09482 + 0.23026 nats.
-    balanced = entropy(torch.tensor([0.5, 0.5]))  # log 2, the largest value on two states.
-    deterministic = entropy(torch.tensor([1.0, 0.0]))  # All mass on one state, so 1 log 1 = 0.
-    uniform_four = entropy(torch.full((4,), 0.25))  # log K at K = 4, proved in Section 1.37.2.
-    show("H((0.9, 0.1))", 0.325, skewed, quiet)
-    show("H((0.5, 0.5)) = log 2", 0.693, balanced, quiet)
-    show("H of a point mass", 0.0, deterministic, quiet)
-    show("H of the uniform on 4 states", 1.386, uniform_four, quiet)
-    return {"skewed": skewed, "balanced": balanced, "uniform_four": uniform_four}
-
-
-# %% 04. Cross-entropy
-def cross_entropies(quiet: bool = False) -> dict:
-    """Section 1.33.3: scoring samples from p with the probabilities of q costs more than H(p)."""
-    cost = cross_entropy(P_REF, Q_MODEL)  # 0.40866 + 0.18326 nats.
-    own = entropy(P_REF)  # 0.17851 + 0.32189 nats, the irreducible part of that cost.
-    same = cross_entropy(P_REF, P_REF)  # H(p, p) recovers the entropy exactly.
-    show("H(p, q)", 0.592, cost, quiet)
-    show("H(p)", 0.500, own, quiet)
-    show("excess cost H(p, q) - H(p)", 0.092, cost - own, quiet)
-    show("H(p, p) - H(p)", 0.0, same - own, quiet)
-    return {"cost": cost, "own": own, "excess": cost - own}
-
-
-# %% 05. Kullback-Leibler divergence
-def kl(samples: int, quiet: bool = False) -> dict:
-    """Section 1.34: the split H(p, q) = H(p) + KL(p || q), nonnegativity, and the asymmetry."""
-    forward = divergence(P_REF, Q_MODEL)  # 0.23015 - 0.13863 nats, the price of using q.
-    reverse = divergence(Q_MODEL, P_REF)  # -0.17261 + 0.27726 nats, the other direction.
-    split = entropy(P_REF) + forward  # The boxed identity, which must return H(p, q).
-    logits = torch.randn(max(samples, 1), 2, 8)  # Random logits for pairs of distributions.
-    pairs = torch.softmax(logits, dim=-1)  # Two distributions over eight states per row; [N, 2, 8].
-    drawn = (pairs[:, 0] * (pairs[:, 0] / pairs[:, 1]).log()).sum(-1)  # KL for every pair; [N].
-    violation = float((-drawn).clamp_min(0.0).max())  # The largest breach of KL >= 0 among them.
-    t = torch.linspace(0.01, 5.0, 500)  # The range on which the inequality log t <= t - 1 is checked.
-    show("KL(p || q)", 0.0915, forward, quiet)
-    show("KL(q || p)", 0.1046, reverse, quiet)
-    show("H(p) + KL(p || q) = H(p, q)", 0.592, split, quiet)
-    show("KL(p || p)", 0.0, divergence(P_REF, P_REF), quiet)
-    show("reverse over forward, percent", 14.0, 100.0 * (reverse / forward - 1.0), quiet)
-    show("largest log t - (t - 1)", 0.0, float((t.log() - (t - 1.0)).max()), quiet)
-    show(f"largest KL breach over {len(drawn)} pairs", 0.0, violation, quiet)
-    return {"forward": forward, "reverse": reverse, "split": split,
-            "smallest": float(drawn.min()), "violation": violation}
-
-
-# %% 06. Cross-entropy by Monte Carlo
-def monte_carlo(samples: int, repeats: int, quiet: bool = False) -> dict:
-    """Section 1.5.3 applied to Section 1.33.3: H(p, q) is an expectation, so an average estimates it.
-
-    We draw states from p and average -log q at the drawn states. The estimate moves with --seed
-    and with --samples, and the exact value it approaches does not. Repeating a ten-sample estimate
-    ``repeats`` times measures the sd that Section 1.5.3 predicts for a ten-sample average.
-    """
-    exact = cross_entropy(P_REF, Q_MODEL)  # The target of the estimator, 0.592 nats.
-    draws = torch.multinomial(P_REF, max(samples, 10), replacement=True)  # [N] states drawn from p.
-    surprise = -Q_MODEL.log()  # The two values -log q(x) the estimator averages over; [2].
-    scored = surprise[draws]  # The summand -log q(x^(n)) at each drawn state; [N].
-    estimate, spread = float(scored.mean()), float(scored.std())  # The average and its per-sample sd.
-    error = spread / math.sqrt(len(draws))  # The sd of the average, which falls like 1/sqrt(N).
-    show("H(p, q) exactly", 0.592, exact, quiet)
-    if not quiet:
-        label = f"H(p, q) from {len(draws)} samples"
-        print(f"  {label:<36s} exact   {exact:>11.5f}   estimate {estimate:>11.5f}"
-              f"   sd {error:.1e}")
-    rounds = max(repeats, 2)  # Two averages is the fewest a standard deviation can be taken over.
-    repeated = torch.multinomial(P_REF, 10 * rounds, replacement=True)  # [10 * rounds] states.
-    averages = surprise[repeated].view(rounds, 10).mean(dim=1)  # One ten-sample estimate per round.
-    exact_spread = float((P_REF @ (surprise - exact).square()).sqrt())  # sd of one scored draw.
-    predicted = exact_spread / math.sqrt(10.0)  # The sd of a ten-sample average, Section 1.5.3.
-    show(f"sd of {rounds} ten-sample estimates", predicted, float(averages.std()), quiet)
-    return {"exact": exact, "estimate": estimate, "error": error, "draws": len(draws),
-            "predicted": predicted, "measured": float(averages.std()), "rounds": rounds}
-
-
-# %% 07. Maximum likelihood as KL minimization
-def maximum_likelihood(quiet: bool = False) -> dict:
-    """Section 1.35: the empirical loss at the frequencies is their entropy, and the gap to any
-    other model is the KL divergence from the frequencies to it."""
-    frequencies = BASE_COUNTS / BASE_COUNTS.sum()  # p_hat = (0.4, 0.3, 0.2, 0.1) over A, C, G, T.
-    uniform = torch.full((4,), 0.25)  # The categorical model that spreads the mass evenly.
-    at_frequencies = cross_entropy(frequencies, frequencies)  # The empirical loss at p_hat.
-    at_uniform = cross_entropy(frequencies, uniform)  # The same data scored by the uniform model.
-    gap = at_uniform - at_frequencies  # The reducible part of the cost, which is a KL divergence.
-    terms = -frequencies * frequencies.log()  # The four terms the chapter adds up; [4].
-    show("term for A, 0.4 log(1/0.4)", 0.36652, float(terms[0]), quiet)
-    show("term for C, 0.3 log(1/0.3)", 0.36119, float(terms[1]), quiet)
-    show("term for G, 0.2 log(1/0.2)", 0.32189, float(terms[2]), quiet)
-    show("term for T, 0.1 log(1/0.1)", 0.23026, float(terms[3]), quiet)
-    show("loss at the frequencies", 1.280, at_frequencies, quiet)
-    show("loss under the uniform model", 1.386, at_uniform, quiet)
-    show("gap between the two losses", 0.1064, gap, quiet)
-    show("KL(p_hat || uniform)", 0.1064, divergence(frequencies, uniform), quiet)
-    return {"frequencies": frequencies, "at_frequencies": at_frequencies,
-            "at_uniform": at_uniform, "gap": gap}
-
-
-# %% 08. The report
-def main() -> None:
-    parser = argparse.ArgumentParser(description=__doc__.split("\n")[0])
-    parser.add_argument("--samples", type=int, default=20000, help="draws in the Monte Carlo checks")
-    parser.add_argument("--steps", type=int, default=5000, help="repeats of the ten-sample estimate")
-    parser.add_argument("--seed", type=int, default=0, help="random seed for the Monte Carlo draws")
-    parser.add_argument("--out", type=str, default="", help="directory for the run's output")
-    parser.add_argument("--quiet", action="store_true", help="print only the final report")
-    args = parser.parse_args()
-
-    set_seed(args.seed)  # Seed first, so the same seed gives the same draws.
-    if not args.quiet:
-        print("Chapter 1, Sections 1.33 to 1.35, worked example by worked example\n")
-    surprise = information(args.quiet)
-    spread = entropies(args.quiet)
-    cost = cross_entropies(args.quiet)
-    gap = kl(args.samples, args.quiet)
-    estimate = monte_carlo(args.samples, args.steps, args.quiet)
-    fitted = maximum_likelihood(args.quiet)
-
-    report = {
-        "I(x) at p(x) = 0.01, nats": surprise["rare"],
-        "H((0.9, 0.1)), nats": spread["skewed"],
-        "H((0.5, 0.5)), nats": spread["balanced"],
-        "H(p), nats": cost["own"],
-        "H(p, q), nats": cost["cost"],
-        "KL(p || q), nats": gap["forward"],
-        "KL(q || p), nats": gap["reverse"],
-        "smallest KL over random pairs": gap["smallest"],
-        "largest KL breach, random pairs": gap["violation"],
-        "H(p, q) from samples": estimate["estimate"],
-        "Monte Carlo samples drawn": estimate["draws"],
-        "ten-sample estimates taken": estimate["rounds"],
-        "sd of those estimates": estimate["measured"],
-        "sd Section 1.5.3 predicts": estimate["predicted"],
-        "empirical loss at p_hat, nats": fitted["at_frequencies"],
-        "loss under the uniform, nats": fitted["at_uniform"],
-        "KL(p_hat || uniform), nats": fitted["gap"],
-    }
-    print_report("Information, Sections 1.33 to 1.35", report)
-    save_run(args.out, config=vars(args), report=report)
-
-
-if __name__ == "__main__":  # Execute when launched directly.
-    main()  # Work every example of Sections 1.33 to 1.35 and print the report.

lecture_1/linear_algebra.py

@@ -1,253 +0,0 @@
-"""The simplex, the softmax, transition matrices and the matrix exponential.
-
-A distribution over K states is a point of the simplex, a softmax maps unrestricted logits onto
-that simplex, and a row-stochastic matrix moves one point of it to the next. Along an
-eigendirection a matrix power is a power of a number, and that is where the stationary
-distribution, the rate at which a chain reaches it, and the matrix exponential of continuous time
-all come from.
-
-Chapter 1, Sections 1.12, 1.13, 1.15, 1.16, 1.17 and 1.18 of the notes. Runs on the three-base
-simplex of Section 1.12.2, the two-state chain of Section 1.15.1 and the symmetric matrix
-A = [[2, 1], [1, 2]] of Section 1.16. Prints the chapter's number beside the computed one.
-
-    python lecture_1/linear_algebra.py             # every worked example
-    python lecture_1/linear_algebra.py --steps 40  # run the chain further toward p*
-    python lecture_1/linear_algebra.py --quiet     # the closing report only
-"""
-
-# %% 01. Imports and configuration
-from __future__ import annotations
-
-import argparse
-import math
-
-import torch
-
-from cis6270.runner import print_report, save_run, set_seed
-
-torch.set_default_dtype(torch.float64)  # Chapter 1 is arithmetic, and five printed digits need it.
-
-
-P_BASES = torch.tensor([0.5, 0.3, 0.2])  # The three-state distribution over A, C, G; Section 1.12.2.
-H_BASES = torch.tensor([1.0, 0.0, 4.0])  # Hydrogen-bond contacts per base, the property h.
-Q_BASES = torch.tensor([0.1, 0.1, 0.8])  # The second point of the simplex, mixed with P_BASES.
-LOGITS = torch.tensor([2.0, 1.0, 0.0])  # The logit vector z of Section 1.13.1.
-P_CHAIN = torch.tensor([[0.8, 0.2],  # Row 1: stay in state one with 0.8, move to state two with 0.2.
-                        [0.3, 0.7]])  # Row 2: move to state one with 0.3, stay with 0.7.
-A_MATRIX = torch.tensor([[2.0, 1.0], [1.0, 2.0]])  # The running matrix of Sections 1.16 and 1.18.
-X_START = torch.tensor([3.0, 1.0])  # The initial vector x_0 = 2 v_1 + v_2 of Section 1.16.2.
-
-
-WORKED: dict[str, float] = {}  # Every computed value, by label, for tests/test_lecture_1.py.
-
-
-def show(label: str, by_hand: float, computed: float, quiet: bool = False) -> None:
-    """Print the chapter's hand-worked value beside the value this file computes, and keep it."""
-    WORKED[label] = computed  # The test suite asserts these against the chapter's own numbers.
-    if not quiet:  # Under --quiet only the closing report is printed.
-        computed = computed + 0.0  # Adding zero prints a negative zero as a plain zero.
-        print(f"  {label:<36s} chapter {by_hand:>11.5f}   computed {computed:>11.5f}"
-              f"   gap {abs(by_hand - computed):.1e}")
-
-
-# %% 02. Probability vectors and the simplex
-def simplex(quiet: bool = False) -> dict:
-    """Sections 1.12.2 and 1.12.3: expectation as a dot product, and a mixture of two points."""
-    expectation = float(P_BASES @ H_BASES)  # E[h(X)] = p^T h = 0.5 + 0 + 0.8 = 1.3.
-    mixture = 0.5 * P_BASES + 0.5 * Q_BASES  # r = lambda p + (1 - lambda) q at lambda = 0.5; [3].
-    corners = 0.3 * torch.tensor([1.0, 0.0]) + 0.7 * torch.tensor([0.0, 1.0])  # Two-state mixture.
-    show("sum_i p_i", 1.0, float(P_BASES.sum()), quiet)
-    show("E[h(X)] = p^T h", 1.3, expectation, quiet)
-    show("mixture r_1 at lambda = 0.5", 0.30, float(mixture[0]), quiet)
-    show("mixture r_2 at lambda = 0.5", 0.20, float(mixture[1]), quiet)
-    show("mixture r_3 at lambda = 0.5", 0.50, float(mixture[2]), quiet)
-    show("sum_i r_i, still a distribution", 1.0, float(mixture.sum()), quiet)
-    show("0.3 e_1 + 0.7 e_2, first entry", 0.3, float(corners[0]), quiet)
-    return {"expectation": expectation, "mixture": mixture}
-
-
-# %% 03. Logits and the softmax
-def softmax(quiet: bool = False) -> dict:
-    """Section 1.13.1: softmax(z)_i = e^{z_i} / sum_j e^{z_j} on the logits z = (2, 1, 0)."""
-    exponentials = LOGITS.exp()  # e^2, e^1, e^0; [3].
-    partition = float(exponentials.sum())  # The normalizing sum, 7.389 + 2.718 + 1.000.
-    probabilities = exponentials / partition  # The point of the simplex the logits name; [3].
-    shifted = torch.softmax(LOGITS + 5.0, dim=0)  # Adding a constant c cancels in the ratio.
-    show("e^{z_1}", 7.389, float(exponentials[0]), quiet)
-    show("e^{z_2}", 2.718, float(exponentials[1]), quiet)
-    show("e^{z_3}", 1.000, float(exponentials[2]), quiet)
-    show("partition function sum_j e^{z_j}", 11.107, partition, quiet)
-    show("softmax(z)_1", 0.665, float(probabilities[0]), quiet)
-    show("softmax(z)_2", 0.245, float(probabilities[1]), quiet)
-    show("softmax(z)_3", 0.090, float(probabilities[2]), quiet)
-    show("sum_i softmax(z)_i", 1.0, float(probabilities.sum()), quiet)
-    show("shift invariance, largest change", 0.0, float((shifted - probabilities).abs().max()), quiet)
-    return {"probabilities": probabilities, "partition": partition}
-
-
-# %% 04. A probability vector under a transition matrix
-def one_transition(quiet: bool = False) -> dict:
-    """Section 1.15.1: p_{t+1} = P^T p_t moves an even split to (0.55, 0.45) in one step."""
-    p_t = torch.tensor([0.5, 0.5])  # The current distribution over the two states; [2].
-    p_next = P_CHAIN.T @ p_t  # Each new entry collects the mass arriving from every state; [2].
-    show("row sums of P, row 1", 1.0, float(P_CHAIN[0].sum()), quiet)
-    show("row sums of P, row 2", 1.0, float(P_CHAIN[1].sum()), quiet)
-    show("p_{t+1}(1)", 0.55, float(p_next[0]), quiet)
-    show("p_{t+1}(2)", 0.45, float(p_next[1]), quiet)
-    show("sum_j p_{t+1}(j)", 1.0, float(p_next.sum()), quiet)
-    return {"p_next": p_next}
-
-
-# %% 05. Eigenvalues and eigenvectors
-def eigen(quiet: bool = False) -> dict:
-    """Section 1.16: the characteristic equation of A, its two eigenvalues and its eigenvectors."""
-    trace = float(A_MATRIX.trace())  # The sum of the eigenvalues, 4.
-    determinant = float(torch.linalg.det(A_MATRIX))  # The product of the eigenvalues, 3.
-    values, vectors = torch.linalg.eigh(A_MATRIX)  # A is symmetric, so the eigenvalues are real.
-    order = values.argsort(descending=True)  # The chapter indexes lambda_1 = 3 before lambda_2 = 1.
-    values, vectors = values[order], vectors[:, order]  # [2] and [2, 2].
-    vectors = vectors / vectors[0]  # Rescale each one to first entry one, which here gives (1, ±1).
-    residual = (A_MATRIX @ vectors - vectors * values).abs().max()  # A v - lambda v at both pairs.
-    show("trace A, lambda_1 + lambda_2", 4.0, trace, quiet)
-    show("det A, lambda_1 lambda_2", 3.0, determinant, quiet)
-    show("lambda_1", 3.0, float(values[0]), quiet)
-    show("lambda_2", 1.0, float(values[1]), quiet)
-    show("v_1 second entry over first", 1.0, float(vectors[1, 0]), quiet)
-    show("v_2 second entry over first", -1.0, float(vectors[1, 1]), quiet)
-    show("largest |A v - lambda v|", 0.0, float(residual), quiet)
-    return {"values": values, "vectors": vectors}
-
-
-# %% 06. Powers along an eigendirection
-def powers(values: torch.Tensor, vectors: torch.Tensor, quiet: bool = False) -> dict:
-    """Section 1.16.2: x_n = c_1 lambda_1^n v_1 + c_2 lambda_2^n v_2, against A^n x_0 taken directly."""
-    coefficients = torch.linalg.solve(vectors, X_START)  # c_1 = 2 and c_2 = 1 for x_0 = (3, 1).
-    eigen_route = vectors @ (coefficients * values ** 2)  # The n = 2 state along the two directions.
-    direct = torch.linalg.matrix_power(A_MATRIX, 2) @ X_START  # A^2 x_0, squaring the matrix first.
-    far = vectors @ (coefficients * values ** 10)  # The same decomposition carried to n = 10.
-    show("c_1", 2.0, float(coefficients[0]), quiet)
-    show("c_2", 1.0, float(coefficients[1]), quiet)
-    show("lambda_1^10", 59049.0, float(values[0]) ** 10, quiet)
-    show("x_2 first entry, eigen route", 19.0, float(eigen_route[0]), quiet)
-    show("x_2 second entry, eigen route", 17.0, float(eigen_route[1]), quiet)
-    show("x_2 first entry, A^2 x_0", 19.0, float(direct[0]), quiet)
-    show("x_10 first entry", 118099.0, float(far[0]), quiet)
-    show("x_10 second entry", 118097.0, float(far[1]), quiet)
-    return {"x_two": eigen_route, "x_ten": far}
-
-
-# %% 07. Diagonalization and the dominant eigenvalue
-def diagonalization(values: torch.Tensor, vectors: torch.Tensor, quiet: bool = False) -> dict:
-    """Sections 1.16.3 and 1.16.4: A = V Lambda V^{-1}, so A^n = V Lambda^n V^{-1} telescopes."""
-    inverse = torch.linalg.inv(vectors)  # V^{-1} = 0.5 [[1, 1], [1, -1]]; [2, 2].
-    rebuilt = vectors @ torch.diag(values) @ inverse  # V Lambda V^{-1} returns A itself; [2, 2].
-    cube = vectors @ torch.diag(values ** 3) @ inverse  # V diag(27, 1) V^{-1}, the n = 3 power.
-    direct_cube = torch.linalg.matrix_power(A_MATRIX, 3)  # A^3 computed by repeated multiplication.
-    dominant = torch.tensor([2.0, 0.5]) ** 10  # The coefficients 2^10 and 0.5^10 of Section 1.16.4.
-    show("V^{-1}_11", 0.5, float(inverse[0, 0]), quiet)
-    show("largest |V Lambda V^{-1} - A|", 0.0, float((rebuilt - A_MATRIX).abs().max()), quiet)
-    show("A^3_11 through the eigenvalues", 14.0, float(cube[0, 0]), quiet)
-    show("A^3_12 through the eigenvalues", 13.0, float(cube[0, 1]), quiet)
-    show("largest |V Lambda^3 V^{-1} - A^3|", 0.0, float((cube - direct_cube).abs().max()), quiet)
-    show("dominant coefficient 2^10", 1024.0, float(dominant[0]), quiet)
-    show("decaying coefficient 0.5^10", 0.000977, float(dominant[1]), quiet)
-    show("ratio of the two, in millions", 1.05, float(dominant[0] / dominant[1]) / 1e6, quiet)
-    return {"inverse": inverse, "cube": cube}
-
-
-# %% 08. Stationary distributions
-def stationary(steps: int, quiet: bool = False) -> dict:
-    """Section 1.17: p* solves (P^T - I) p* = 0 with sum one, and the chain halves its deviation."""
-    drift = (P_CHAIN.T - torch.eye(2))[0]  # The first row of P^T - I, which reads 0.3 p_2 = 0.2 p_1.
-    system = torch.stack([drift, torch.ones(2)])  # That one equation, with normalization; [2, 2].
-    p_star = torch.linalg.solve(system, torch.tensor([0.0, 1.0]))  # The lambda = 1 eigenvector; [2].
-    trace, determinant = float(P_CHAIN.trace()), float(torch.linalg.det(P_CHAIN))  # 1.5 and 0.5.
-    second = (trace - math.sqrt(trace ** 2 - 4.0 * determinant)) / 2.0  # The smaller root, 0.5.
-    p_t, history = torch.tensor([0.5, 0.5]), []  # Start from the even split of Section 1.15.1.
-    for _ in range(max(steps, 1)):  # Iterate p_{t+1} = P^T p_t for as many steps as requested.
-        p_t = P_CHAIN.T @ p_t  # One more application of the transition matrix; [2].
-        history.append(p_t.clone())
-    show("p*_1", 0.6, float(p_star[0]), quiet)
-    show("p*_2", 0.4, float(p_star[1]), quiet)
-    show("largest |P^T p* - p*|", 0.0, float((P_CHAIN.T @ p_star - p_star).abs().max()), quiet)
-    show("trace P", 1.5, trace, quiet)
-    show("det P", 0.5, determinant, quiet)
-    show("second eigenvalue, the rate", 0.5, second, quiet)
-    for index, value in ((0, 0.55), (1, 0.575), (2, 0.5875), (9, 0.59990)):
-        if index < len(history):  # The chapter prints p_1, p_2, p_3 and p_10 of this sequence.
-            show(f"p_{index + 1}(1)", value, float(history[index][0]), quiet)
-    deviation = float((history[-1] - p_star).abs().max())  # |d_n| = 0.5^n |d_0| after n steps.
-    show(f"|d_n| after {len(history)} steps", 0.1 * second ** len(history), deviation, quiet)
-    return {"p_star": p_star, "second": second, "history": history, "deviation": deviation}
-
-
-# %% 09. From discrete steps to continuous time
-def matrix_exponential(values: torch.Tensor, vectors: torch.Tensor, quiet: bool = False) -> dict:
-    """Section 1.18.3: the four-term series and the eigenvalue formula for e^{tA} at t = 0.1."""
-    t = 0.1  # The small time at which the chapter truncates the series.
-    series, power, factorial = torch.eye(2), torch.eye(2), 1.0  # Accumulate I + tA + t^2A^2/2 + ...
-    for k in range(1, 4):  # Terms k = 1, 2, 3, which with the identity make four.
-        power, factorial = power @ A_MATRIX, factorial * k  # A^k and k factorial.
-        series = series + t ** k * power / factorial  # Add t^k A^k / k!; [2, 2].
-    discarded = t ** 4 * (power @ A_MATRIX) / 24.0  # The first omitted term, t^4 A^4 / 4!.
-    scaled = torch.diag((t * values).exp())  # e^{t Lambda}, the diagonal of entry exponentials; [2, 2].
-    eigen_route = vectors @ scaled @ torch.linalg.inv(vectors)  # V e^{t Lambda} V^{-1}; [2, 2].
-    gap = float((series - eigen_route).abs().max())  # How far four terms fall short at t = 0.1.
-    coefficients = torch.linalg.solve(vectors, X_START)  # c_1 = 2, c_2 = 1, as in Section 1.16.2.
-    x_one = vectors @ (coefficients * values.exp())  # x_t = sum_i c_i e^{t lambda_i} v_i at t = 1.
-    show("four-term series, entry 11", 1.22733, float(series[0, 0]), quiet)
-    show("four-term series, entry 12", 0.12217, float(series[0, 1]), quiet)
-    show("eigenvalue route, entry 11", 1.22751, float(eigen_route[0, 0]), quiet)
-    show("eigenvalue route, entry 12", 0.12234, float(eigen_route[0, 1]), quiet)
-    show("largest gap between the routes", 1.8e-4, gap, quiet)
-    show("first discarded term, entry 11", 1.71e-4, float(discarded[0, 0]), quiet)
-    show("first discarded term, entry 12", 1.67e-4, float(discarded[0, 1]), quiet)
-    show("x_1 first entry, 2e^3 + e", 42.889, float(x_one[0]), quiet)
-    show("x_1 second entry, 2e^3 - e", 37.453, float(x_one[1]), quiet)
-    show("decay e^{-2t} at t = 1", 0.1353, math.exp(-2.0), quiet)
-    show("decay e^{-2t} at t = 3", 0.002479, math.exp(-6.0), quiet)
-    show("sqrt(lambda_1), Gaussian long axis", 1.732, math.sqrt(float(values[0])), quiet)
-    return {"series": series, "eigen_route": eigen_route, "gap": gap, "x_one": x_one}
-
-
-# %% 10. The report
-def main() -> None:
-    parser = argparse.ArgumentParser(description=__doc__.split("\n")[0])
-    parser.add_argument("--steps", type=int, default=10, help="iterations of p_{t+1} = P^T p_t")
-    parser.add_argument("--seed", type=int, default=0, help="random seed, for a repeatable run")
-    parser.add_argument("--out", type=str, default="", help="directory for the run's output")
-    parser.add_argument("--quiet", action="store_true", help="print only the final report")
-    args = parser.parse_args()
-
-    set_seed(args.seed)  # Seed first, so the same seed gives the same run.
-    if not args.quiet:
-        print("Chapter 1, Sections 1.12 to 1.18, worked example by worked example\n")
-    triangle = simplex(args.quiet)
-    logits = softmax(args.quiet)
-    step = one_transition(args.quiet)
-    spectrum = eigen(args.quiet)
-    grown = powers(spectrum["values"], spectrum["vectors"], args.quiet)
-    diagonalization(spectrum["values"], spectrum["vectors"], args.quiet)
-    chain = stationary(args.steps, args.quiet)
-    flow = matrix_exponential(spectrum["values"], spectrum["vectors"], args.quiet)
-
-    report = {
-        "E[h(X)] on the three bases": triangle["expectation"],
-        "softmax(z)_1": float(logits["probabilities"][0]),
-        "p_{t+1}(1) after one step": float(step["p_next"][0]),
-        "lambda_1 of A": float(spectrum["values"][0]),
-        "lambda_2 of A": float(spectrum["values"][1]),
-        "x_10 first entry": float(grown["x_ten"][0]),
-        "p*_1, stationary": float(chain["p_star"][0]),
-        "second eigenvalue of P": chain["second"],
-        "chain steps run": len(chain["history"]),
-        "deviation from p* at the end": chain["deviation"],
-        "e^{tA}_11 at t = 0.1": float(flow["eigen_route"][0, 0]),
-        "four-term truncation error": flow["gap"],
-    }
-    print_report("Linear algebra, Sections 1.12 to 1.18", report)
-    save_run(args.out, config=vars(args), report=report)
-
-
-if __name__ == "__main__":  # Execute when launched directly.
-    main()  # Work every example of Sections 1.12 to 1.18 and print the report.

lecture_1/optimization.py

@@ -1,223 +0,0 @@
-"""Unconstrained descent, Lagrange multipliers, maximum entropy and the KKT conditions.
-
-A differentiable minimum has a zero gradient, and a descent step moves against that gradient. An
-equality constraint comes with a multiplier, and at the constrained optimum the two gradients are
-parallel. Applied to entropy under an expectation constraint, that condition leads to an
-exponential family. An inequality comes with a sign requirement and complementary slackness, and
-under the softmax the simplex constraints hold by construction.
-
-Chapter 1, Sections 1.36, 1.37 and 1.38 of the notes. Runs on the quadratic loss of
-Section 1.36, the constrained distance of Section 1.37.1, the GC-content target of
-Section 1.37.3 and the bounded quadratic of Section 1.38.
-
-    python lecture_1/optimization.py              # every worked example
-    python lecture_1/optimization.py --steps 20   # a short descent, far from the minimizer
-    python lecture_1/optimization.py --quiet      # the closing report only
-"""
-
-# %% 01. Imports and configuration
-from __future__ import annotations
-
-import argparse
-import math
-
-import torch
-
-from cis6270.runner import print_report, save_run, set_seed
-
-torch.set_default_dtype(torch.float64)  # Chapter 1 is arithmetic, and five printed digits need it.
-
-
-BASE_COUNTS = torch.tensor([4.0, 3.0, 2.0, 1.0])  # n(A), n(C), n(G), n(T) over ten DNA bases.
-GC_INDICATOR = torch.tensor([0.0, 1.0, 1.0, 0.0])  # h(x) = 1 on the strong bases C and G.
-GC_TARGET = 0.6  # The required expected GC content, the constraint value c of Section 1.37.3.
-
-
-WORKED: dict[str, float] = {}  # Every computed value, by label, for tests/test_lecture_1.py.
-
-
-def show(label: str, by_hand: float, computed: float, quiet: bool = False) -> None:
-    """Print the chapter's hand-worked value beside the value this file computes, and keep it."""
-    WORKED[label] = computed  # The test suite asserts these against the chapter's own numbers.
-    if not quiet:  # Under --quiet only the closing report is printed.
-        computed = computed + 0.0  # Adding zero prints a negative zero as a plain zero.
-        print(f"  {label:<36s} chapter {by_hand:>11.5f}   computed {computed:>11.5f}"
-              f"   gap {abs(by_hand - computed):.1e}")
-
-
-def entropy(p: torch.Tensor) -> float:
-    """H(p) = -sum_i p_i log p_i, the quantity maximized in both constrained problems below."""
-    return float(-(p * p.clamp_min(1e-300).log()).sum())  # Clamp keeps 0 log 0 at zero.
-
-
-# %% 02. Optimization without constraints
-def unconstrained(steps: int, quiet: bool = False) -> dict:
-    """Section 1.36: on L(theta) = (theta - 3)^2 the error contracts by 0.8 at every step."""
-    learning_rate, theta = 0.1, 0.0  # The step size eta and the starting parameter theta_0.
-    history = []  # The iterates, so the first three can be compared with the chapter.
-    for _ in range(max(steps, 1)):
-        theta = theta - learning_rate * 2.0 * (theta - 3.0)  # theta - eta L'(theta) = 0.8 theta + 0.6.
-        history.append(theta)
-    closed_form = 3.0 - 3.0 * 0.8 ** len(history)  # theta_k = 3 - 3(0.8)^k, from the contraction.
-    for index, value in ((0, 0.600), (1, 1.080), (2, 1.464)):
-        if index < len(history):  # The chapter prints theta_1, theta_2 and theta_3.
-            show(f"theta_{index + 1}", value, history[index], quiet)
-    show(f"theta_{len(history)} from the closed form", closed_form, history[-1], quiet)
-    show("gradient at the minimizer theta = 3", 0.0, 2.0 * (3.0 - 3.0), quiet)
-    return {"history": history, "closed_form": closed_form}
-
-
-# %% 03. An equality constraint and its multiplier
-def equality_constraint(quiet: bool = False) -> dict:
-    """Section 1.37.1: minimize x_1^2 + x_2^2 on the line x_1 + x_2 = 1, through the Lagrangian."""
-    system = torch.tensor([[2.0, 0.0, 1.0],  # 2 x_1 + lambda = 0, the first stationarity condition.
-                           [0.0, 2.0, 1.0],  # 2 x_2 + lambda = 0, the second.
-                           [1.0, 1.0, 0.0]])  # x_1 + x_2 = 1, the constraint itself.
-    solution = torch.linalg.solve(system, torch.tensor([0.0, 0.0, 1.0]))  # (x_1, x_2, lambda); [3].
-    optimum = float(solution[:2].square().sum())  # f at the constrained minimizer.
-    endpoint = float(torch.tensor([1.0, 0.0]).square().sum())  # f at either end of the line.
-    grad_f = 2.0 * solution[:2]  # grad f = (2 x_1, 2 x_2) at the solution; [2].
-    parallel = float((grad_f + solution[2] * torch.ones(2)).abs().max())  # grad f + lambda grad g.
-    show("x_1 at the constrained minimum", 0.5, float(solution[0]), quiet)
-    show("x_2 at the constrained minimum", 0.5, float(solution[1]), quiet)
-    show("lambda, the equality multiplier", -1.0, float(solution[2]), quiet)
-    show("f at the minimum", 0.5, optimum, quiet)
-    show("f at the endpoint (1, 0)", 1.0, endpoint, quiet)
-    show("largest |grad f + lambda grad g|", 0.0, parallel, quiet)
-    return {"solution": solution, "optimum": optimum}
-
-
-# %% 04. Normalization as a constraint
-def normalization(quiet: bool = False) -> dict:
-    """Section 1.37.2: entropy under normalization alone is maximized by the uniform distribution,
-    and the empirical negative log-likelihood is minimized by the observed frequencies."""
-    states = 4  # K, the size of the DNA alphabet.
-    uniform = torch.full((states,), 1.0 / states)  # p_i = e^{lambda - 1}, the same at every state.
-    frequencies = BASE_COUNTS / BASE_COUNTS.sum()  # p(x) = n(x) / N, the solution at lambda = 1.
-    loss = float(-(frequencies * frequencies.log()).sum())  # The minimum empirical loss, its entropy.
-    nearby = frequencies + torch.tensor([0.05, -0.05, 0.0, 0.0])  # A different categorical model.
-    worse = float(-(frequencies * nearby.log()).sum())  # Its loss on the same ten bases.
-    show("uniform p_i at K = 4", 0.25, float(uniform[0]), quiet)
-    show("max entropy under normalization", 1.386, entropy(uniform), quiet)
-    show("p_hat(A) = n(A) / N", 0.4, float(frequencies[0]), quiet)
-    show("p_hat(T) = n(T) / N", 0.1, float(frequencies[3]), quiet)
-    show("minimum empirical loss", 1.280, loss, quiet)
-    show("loss drop from moving off p_hat", 0.0, max(0.0, loss - worse), quiet)
-    if not quiet:  # The displaced model is not a chapter number, so it is printed on its own.
-        print(f"  {'loss of a model moved by 0.05':<36s} loss    {worse:>11.5f}"
-              f"   at p_hat {loss:>11.5f}")
-    return {"uniform": uniform, "frequencies": frequencies, "loss": loss, "worse": worse}
-
-
-# %% 05. Maximum entropy under a property constraint
-def maximum_entropy(samples: int, quiet: bool = False) -> dict:
-    """Section 1.37.3: with one multiplier per property the entropy maximizer is p* proportional to
-    exp(lambda h(x)), and the constraint fixes lambda."""
-    strong, weak = float(GC_INDICATOR.sum()), float((1.0 - GC_INDICATOR).sum())  # Two bases each.
-    exp_lambda = (GC_TARGET / (1.0 - GC_TARGET)) * (weak / strong)  # From 2e^L/(2e^L + 2) = 0.6.
-    multiplier = math.log(exp_lambda)  # lambda = log 1.5, the tilt toward the strong bases.
-    unnormalized = (multiplier * GC_INDICATOR).exp()  # exp(lambda h(x)) at the four bases; [4].
-    partition = float(unnormalized.sum())  # Z(lambda) = 2 e^lambda + 2 = 5.
-    p_star = unnormalized / partition  # The maximum-entropy distribution (0.2, 0.3, 0.3, 0.2); [4].
-    uniform = torch.full((4,), 0.25)  # What lambda = 0 returns, with no property requirement.
-    weights = torch.rand(max(samples, 1), 2)  # One split per base class, drawn independently; [N, 2].
-    competitors = torch.stack([  # Distributions that also meet the GC target, in A, C, G, T order.
-        (1.0 - GC_TARGET) * weights[:, 0], GC_TARGET * weights[:, 1],
-        GC_TARGET * (1.0 - weights[:, 1]), (1.0 - GC_TARGET) * (1.0 - weights[:, 0])], dim=1)
-    drawn = -(competitors * competitors.clamp_min(1e-300).log()).sum(1)  # H of every competitor; [N].
-    best = float(drawn.max())  # The largest entropy reached by any distribution meeting the target.
-    show("e^lambda", 1.5, exp_lambda, quiet)
-    show("lambda = log 1.5", 0.4055, multiplier, quiet)
-    show("partition function Z", 5.0, partition, quiet)
-    show("p*(A)", 0.2, float(p_star[0]), quiet)
-    show("p*(C)", 0.3, float(p_star[1]), quiet)
-    show("GC content of p*", 0.6, float(p_star @ GC_INDICATOR), quiet)
-    show("H(p*)", 1.366, entropy(p_star), quiet)
-    show("H of the unconstrained uniform", 1.386, entropy(uniform), quiet)
-    show("entropy cost of the constraint", 0.0201, entropy(uniform) - entropy(p_star), quiet)
-    show(f"H above H(p*) over {len(drawn)} draws", 0.0, max(0.0, best - entropy(p_star)), quiet)
-    return {"multiplier": multiplier, "p_star": p_star, "partition": partition, "best_drawn": best}
-
-
-# %% 06. Inequality constraints and the KKT conditions
-def kkt(quiet: bool = False) -> dict:
-    """Section 1.38: on min (x - 3)^2 with x <= 1 the bound is active, so lambda is positive."""
-    x_star = 1.0  # The unconstrained minimizer 3 breaks the bound, so the solution is at it.
-    multiplier = -2.0 * (x_star - 3.0)  # Stationarity 2(x - 3) + lambda = 0 solved at x = 1.
-    slack = x_star - 1.0  # g(x*) = x* - 1, which must be at or below zero.
-    show("x* at the bound", 1.0, x_star, quiet)
-    show("lambda, the inequality multiplier", 4.0, multiplier, quiet)
-    show("g(x*) = x* - 1", 0.0, slack, quiet)
-    show("complementary slackness lambda g", 0.0, multiplier * slack, quiet)
-    show("stationarity 2(x* - 3) + lambda", 0.0, 2.0 * (x_star - 3.0) + multiplier, quiet)
-    show("f at the unconstrained minimum", 0.0, (3.0 - 3.0) ** 2, quiet)
-    show("f under the bound", 4.0, (x_star - 3.0) ** 2, quiet)
-    show("value drop for epsilon = 0.001", 0.004, -((1.001 - 3.0) ** 2 - 4.0), quiet)
-    return {"multiplier": multiplier, "value": (x_star - 3.0) ** 2}
-
-
-# %% 07. Optimization over the simplex through the softmax
-def over_the_simplex(steps: int, quiet: bool = False) -> dict:
-    """Section 1.38.1: under logits the simplex constraints hold at every step, so the constrained
-    likelihood problem of Section 1.37.2 becomes an unconstrained descent that reaches p_hat."""
-    frequencies = BASE_COUNTS / BASE_COUNTS.sum()  # The minimizer the multiplier route produced.
-    logits = torch.zeros(4, requires_grad=True)  # z in R^4, starting at the uniform distribution.
-    learning_rate, losses = 1.0, []  # Plain gradient descent on z, with the loss after each step.
-    for _ in range(max(steps, 1)):
-        model = torch.softmax(logits, dim=0)  # p_i = e^{z_i} / sum_j e^{z_j}, always on the simplex.
-        loss = -(frequencies * model.log()).sum()  # The empirical negative log-likelihood.
-        gradient, = torch.autograd.grad(loss, logits)  # d L / d z, which is softmax(z) - p_hat.
-        with torch.no_grad():  # The update itself needs no graph.
-            logits -= learning_rate * gradient  # One descent step of Section 1.36, now on logits.
-        losses.append(float(loss.detach()))
-    final = torch.softmax(logits.detach(), dim=0)  # The distribution the descent arrived at; [4].
-    show("sum_i p_i at the end", 1.0, float(final.sum()), quiet)
-    show("entries at or below zero", 0.0, float((final <= 0.0).sum()), quiet)
-    show("p(A) after the descent", 0.4, float(final[0]), quiet)
-    show("p(T) after the descent", 0.1, float(final[3]), quiet)
-    show("loss after the descent", 1.280, losses[-1], quiet)
-    show("largest |p - p_hat|", 0.0, float((final - frequencies).abs().max()), quiet)
-    return {"final": final, "losses": losses, "distance": float((final - frequencies).abs().max())}
-
-
-# %% 08. The report
-def main() -> None:
-    parser = argparse.ArgumentParser(description=__doc__.split("\n")[0])
-    parser.add_argument("--steps", type=int, default=600, help="gradient descent steps")
-    parser.add_argument("--samples", type=int, default=4000, help="competing distributions drawn")
-    parser.add_argument("--seed", type=int, default=0, help="random seed for the drawn competitors")
-    parser.add_argument("--out", type=str, default="", help="directory for the run's output")
-    parser.add_argument("--quiet", action="store_true", help="print only the final report")
-    args = parser.parse_args()
-
-    set_seed(args.seed)  # Seed first, so the same seed gives the same competitors.
-    if not args.quiet:
-        print("Chapter 1, Sections 1.36 to 1.38, worked example by worked example\n")
-    descent = unconstrained(args.steps, args.quiet)
-    constrained = equality_constraint(args.quiet)
-    normalized = normalization(args.quiet)
-    tilted = maximum_entropy(args.samples, args.quiet)
-    bounded = kkt(args.quiet)
-    fitted = over_the_simplex(args.steps, args.quiet)
-
-    report = {
-        "descent steps run": len(descent["history"]),
-        "theta after the descent": descent["history"][-1],
-        "x_1 under x_1 + x_2 = 1": float(constrained["solution"][0]),
-        "lambda, equality constraint": float(constrained["solution"][2]),
-        "minimum empirical loss, nats": normalized["loss"],
-        "lambda for GC content 0.6": tilted["multiplier"],
-        "p*(C) under that constraint": float(tilted["p_star"][1]),
-        "H(p*), nats": entropy(tilted["p_star"]),
-        "best drawn H at the same target": tilted["best_drawn"],
-        "lambda, active inequality": bounded["multiplier"],
-        "f under the bound": bounded["value"],
-        "loss after the softmax descent": fitted["losses"][-1],
-        "largest |p - p_hat|": fitted["distance"],
-    }
-    print_report("Optimization, Sections 1.36 to 1.38", report)
-    save_run(args.out, config=vars(args), report=report)
-
-
-if __name__ == "__main__":  # Execute when launched directly.
-    main()  # Work every example of Sections 1.36 to 1.38 and print the report.

lecture_1/probability.py

@@ -1,261 +0,0 @@
-"""Probability on discrete and continuous states, in the chapter's own numbers.
-
-A distribution allocates one unit of mass over the states, expectation averages a property under
-that allocation, variance and covariance record how the values fluctuate, and conditioning and
-Bayes' rule revise the allocation once something has been observed. The chapter works every
-example below by hand, so each section prints the hand value beside the value this file computes.
-
-Chapter 1, Sections 1.3 to 1.11 of the notes. Runs on the two-token DNA distribution of
-Section 1.3.1, the three peptides of Section 1.7.1 and the motif-and-binding table of
-Section 1.9.1. Prints the chapter's number, the computed number and the gap between them.
-
-    python lecture_1/probability.py                   # every worked example
-    python lecture_1/probability.py --samples 100000  # a longer Monte Carlo ladder
-    python lecture_1/probability.py --quiet           # the closing report only
-"""
-
-# %% 01. Imports and configuration
-from __future__ import annotations
-
-import argparse
-import math
-
-import torch
-
-from cis6270.runner import print_report, save_run, set_seed
-
-torch.set_default_dtype(torch.float64)  # Chapter 1 is arithmetic, and five printed digits need it.
-
-
-STATES = ["AA", "AC", "CA", "CC"]  # The state space X of Section 1.3.1, with L = 2 over {A, C}.
-P_X = torch.tensor([0.4, 0.1, 0.2, 0.3])  # p_X over those four states, in the order of STATES.
-A_COUNT = torch.tensor([2.0, 1.0, 1.0, 0.0])  # h(x), the number of A tokens in each state.
-PEPTIDES = torch.tensor([[1.0, 2.0], [2.0, 3.0], [3.0, 7.0]])  # (hydrophobicity, binding score).
-JOINT = torch.tensor([[0.32, 0.08],  # Row x = 0, motif absent: columns are y = 0 and y = 1.
-                      [0.24, 0.36]])  # Row x = 1, motif present: Table 1.1 of Section 1.9.1.
-
-
-WORKED: dict[str, float] = {}  # Every computed value, by label, for tests/test_lecture_1.py.
-
-
-def show(label: str, by_hand: float, computed: float, quiet: bool = False) -> None:
-    """Print the chapter's hand-worked value beside the value this file computes, and keep it."""
-    WORKED[label] = computed  # The test suite asserts these against the chapter's own numbers.
-    if not quiet:  # Under --quiet only the closing report is printed.
-        computed = computed + 0.0  # Adding zero prints a negative zero as a plain zero.
-        print(f"  {label:<36s} chapter {by_hand:>11.5f}   computed {computed:>11.5f}"
-              f"   gap {abs(by_hand - computed):.1e}")
-
-
-# %% 02. A distribution over length-two sequences
-def two_token_distribution(quiet: bool = False) -> dict:
-    """Section 1.3.1: the four masses sum to one, and the event {x_1 = A} has probability 0.5."""
-    total = float(P_X.sum())  # Normalization, sum over x in X of p_X(x).
-    event = [index for index, state in enumerate(STATES) if state[0] == "A"]  # A = {x : x_1 = A}.
-    first_is_a = float(P_X[event].sum())  # P(X in A), the masses of {AA, AC} added together.
-    show("sum_x p_X(x)", 1.0, total, quiet)
-    show("P(x_1 = A) = p(AA) + p(AC)", 0.5, first_is_a, quiet)
-    return {"normalization": total, "p_first_is_a": first_is_a}
-
-
-# %% 03. Densities for continuous states
-def continuous_densities(quiet: bool = False) -> dict:
-    """Section 1.4: a density is probability per unit length, so its height may exceed one."""
-    uniform_height = 1.0 / 0.1  # The uniform density on [0, 0.1] has the value 1/0.1 there.
-    narrow_mass = uniform_height * 0.02  # P(0 <= X <= 0.02) is the height times the width.
-    gaussian_peak = 1.0 / (math.sqrt(2.0 * math.pi) * 0.1)  # N(0, 0.1^2) at its center.
-    show("uniform[0, 0.1] density height", 10.0, uniform_height, quiet)
-    show("P(0 <= X <= 0.02)", 0.2, narrow_mass, quiet)
-    show("peak of N(0, 0.1^2)", 3.989, gaussian_peak, quiet)
-    return {"uniform_height": uniform_height, "gaussian_peak": gaussian_peak}
-
-
-# %% 04. The expected adenine count
-def expected_adenine(quiet: bool = False) -> dict:
-    """Section 1.5.1: E[h(X)] = sum_x p_X(x) h(x), the probability-weighted average of the count."""
-    expectation = float(P_X @ A_COUNT)  # 0.4(2) + 0.1(1) + 0.2(1) + 0.3(0).
-    across_thousand = 1000.0 * expectation  # The adenines expected in 1000 generated sequences.
-    show("E[h(X)], expected A count", 1.1, expectation, quiet)
-    show("expected A count in 1000 draws", 1100.0, across_thousand, quiet)
-    return {"expectation": expectation}
-
-
-# %% 05. Monte Carlo estimation
-def monte_carlo(samples: int, repeats: int, quiet: bool = False) -> dict:
-    """Section 1.5.3: the sample average estimates E[h(X)] with an error that shrinks like 1/sqrt(N).
-
-    The chapter's own ten-sample draw gave 1.2 against the exact 1.1, a gap well inside the
-    standard deviation 0.263 of a ten-sample average. We redraw here, so the estimate moves with
-    --seed while the exact value and the error scale stay put. Repeating the ten-sample estimate
-    ``repeats`` times turns that predicted 0.263 into a measured spread.
-    """
-    exact = float(P_X @ A_COUNT)  # The target of the estimator, E[h(X)] = 1.1.
-    second_moment = float(P_X @ A_COUNT.square())  # E[h^2] = 0.4(4) + 0.1 + 0.2 = 1.9.
-    variance = second_moment - exact ** 2  # Var(h(X)) = 1.9 - 1.21 = 0.69.
-    show("E[h(X)^2]", 1.9, second_moment, quiet)
-    show("Var(h(X))", 0.69, variance, quiet)
-    ladder, size = [], 10  # The chapter quotes the error at N = 10 and at N = 1000.
-    while size <= max(samples, 10):
-        draws = torch.multinomial(P_X, size, replacement=True)  # [size] state indices from p_X.
-        estimate = float(A_COUNT[draws].mean())  # (1/N) sum_n h(x^(n)), the Monte Carlo average.
-        ladder.append((size, estimate, math.sqrt(variance / size)))
-        size *= 10
-    if not quiet:
-        for size, estimate, error in ladder:
-            print(f"  N = {size:<8d} estimate {estimate:.4f}   exact {exact:.4f}"
-                  f"   sd of the average {error:.4f}")
-    show("sd of a 10-sample average", 0.263, math.sqrt(variance / 10), quiet)
-    show("sd of a 1000-sample average", 0.0263, math.sqrt(variance / 1000), quiet)
-    rounds = max(repeats, 2)  # Two averages is the fewest a standard deviation can be taken over.
-    repeated = torch.multinomial(P_X, 10 * rounds, replacement=True)  # [10 * rounds] state indices.
-    averages = A_COUNT[repeated].view(rounds, 10).mean(dim=1)  # One ten-sample average per round.
-    measured = float(averages.std())  # The spread of those averages, which 0.263 predicts.
-    show(f"sd of {rounds} ten-sample averages", 0.263, measured, quiet)
-    return {"exact": exact, "variance": variance, "ladder": ladder, "measured": measured,
-            "rounds": rounds}
-
-
-# %% 06. A worked variance
-def worked_variance(quiet: bool = False) -> dict:
-    """Section 1.6.2: Var(X) = E[X^2] - (E[X])^2 on the values 0, 1, 2 with masses 0.2, 0.5, 0.3."""
-    values = torch.tensor([0.0, 1.0, 2.0])  # The three values X can take.
-    masses = torch.tensor([0.2, 0.5, 0.3])  # p_X at those three values.
-    first_moment = float(masses @ values)  # E[X] = 0 + 0.5 + 0.6 = 1.1.
-    second_moment = float(masses @ values.square())  # E[X^2] = 0 + 0.5 + 1.2 = 1.7.
-    from_moments = second_moment - first_moment ** 2  # The boxed identity of Section 1.6.1.
-    from_definition = float(masses @ (values - first_moment).square())  # The two-pass route.
-    show("E[X]", 1.1, first_moment, quiet)
-    show("E[X^2]", 1.7, second_moment, quiet)
-    show("Var(X) from the two moments", 0.49, from_moments, quiet)
-    show("Var(X) from the definition", 0.49, from_definition, quiet)
-    show("sigma_X", 0.7, math.sqrt(from_moments), quiet)
-    return {"variance": from_moments, "sigma": math.sqrt(from_moments)}
-
-
-# %% 07. Covariance, correlation and the covariance matrix
-def peptide_covariance(quiet: bool = False) -> dict:
-    """Sections 1.7.1 and 1.7.2: three equally likely peptides, their correlation and their Sigma."""
-    mu = PEPTIDES.mean(dim=0)  # (mu_X, mu_Y) = (2, 4), each coordinate mean; [2].
-    centered = PEPTIDES - mu  # The deviations (-1, -2), (0, -1), (1, 3); [3, 2].
-    sigma_matrix = centered.T @ centered / PEPTIDES.shape[0]  # E[(X-mu)(X-mu)^T]; [2, 2].
-    cov_xy = float(sigma_matrix[0, 1])  # Cov(X, Y) = (2 + 0 + 3)/3 = 5/3.
-    sd_x = math.sqrt(float(sigma_matrix[0, 0]))  # sigma_X = sqrt(2/3).
-    sd_y = math.sqrt(float(sigma_matrix[1, 1]))  # sigma_Y = sqrt(14/3).
-    rho = cov_xy / (sd_x * sd_y)  # Correlation, which divides out the units of both variables.
-    show("mu_X", 2.0, float(mu[0]), quiet)
-    show("mu_Y", 4.0, float(mu[1]), quiet)
-    show("Cov(X, Y)", 5.0 / 3.0, cov_xy, quiet)
-    show("Var(X) = Sigma_11", 2.0 / 3.0, float(sigma_matrix[0, 0]), quiet)
-    show("Var(Y) = Sigma_22", 14.0 / 3.0, float(sigma_matrix[1, 1]), quiet)
-    show("sigma_X across the peptides", 0.8165, sd_x, quiet)
-    show("sigma_Y across the peptides", 2.160, sd_y, quiet)
-    show("rho_XY", 0.9449, rho, quiet)
-    show("Sigma symmetry, Sigma_21 - Sigma_12", 0.0, float(sigma_matrix[1, 0] - cov_xy), quiet)
-    return {"sigma_matrix": sigma_matrix, "rho": rho}
-
-
-# %% 08. Gaussian noise and the noising equation
-def noising_covariance(quiet: bool = False) -> dict:
-    """Section 1.8: Cov(X_t) = alpha_t^2 Cov(X_0) + sigma_t^2 I at alpha_t = 0.6, sigma_t = 0.8."""
-    alpha_t, sigma_t = 0.6, 0.8  # A variance-preserving pair, since 0.36 + 0.64 = 1.
-    identity = torch.eye(2)  # The isotropic noise covariance; [2, 2].
-    isotropic = alpha_t ** 2 * identity + sigma_t ** 2 * identity  # Data covariance I stays I.
-    wide_data = torch.diag(torch.tensor([4.0, 1.0]))  # Cov(X_0) = diag(4, 1); [2, 2].
-    anisotropic = alpha_t ** 2 * wide_data + sigma_t ** 2 * identity  # diag(2.08, 1.00).
-    show("alpha_t^2 + sigma_t^2", 1.0, alpha_t ** 2 + sigma_t ** 2, quiet)
-    show("Cov(X_t)_11 with Cov(X_0) = I", 1.0, float(isotropic[0, 0]), quiet)
-    show("Cov(X_t)_11 with diag(4, 1)", 2.08, float(anisotropic[0, 0]), quiet)
-    show("Cov(X_t)_22 with diag(4, 1)", 1.00, float(anisotropic[1, 1]), quiet)
-    return {"anisotropic": anisotropic}
-
-
-# %% 09. Joint, marginal and conditional distributions
-def joint_marginal_conditional(quiet: bool = False) -> dict:
-    """Sections 1.9.1 to 1.9.3: the motif-and-binding table, its marginals and its conditionals."""
-    total = float(JOINT.sum())  # 0.36 + 0.24 + 0.08 + 0.32 = 1.
-    p_x = JOINT.sum(dim=1)  # p(x) = sum_y p(x, y), motif present or absent; [2].
-    p_y = JOINT.sum(dim=0)  # p(y) = sum_x p(x, y), high or low binding; [2].
-    p_y_given_x = JOINT / p_x[:, None]  # p(y | x) = p(x, y) / p(x), renormalized row by row; [2, 2].
-    show("sum over the four cells", 1.0, total, quiet)
-    show("p(X = 1), motif present", 0.60, float(p_x[1]), quiet)
-    show("p(X = 0), motif absent", 0.40, float(p_x[0]), quiet)
-    show("p(Y = 1), high binding", 0.44, float(p_y[1]), quiet)
-    show("p(Y = 0), low binding", 0.56, float(p_y[0]), quiet)
-    show("p(Y = 1 | X = 1)", 0.6, float(p_y_given_x[1, 1]), quiet)
-    show("p(Y = 0 | X = 1)", 0.4, float(p_y_given_x[1, 0]), quiet)
-    show("p(Y = 1 | X = 0)", 0.2, float(p_y_given_x[0, 1]), quiet)
-    show("p(Y = 0 | X = 0)", 0.8, float(p_y_given_x[0, 0]), quiet)
-    return {"p_x": p_x, "p_y": p_y, "p_y_given_x": p_y_given_x}
-
-
-# %% 10. Bayes' rule and the law of total probability
-def bayes(quiet: bool = False) -> dict:
-    """Section 1.10: the motif probability after observing high binding, with its own denominator."""
-    p_x = JOINT.sum(dim=1)  # The prior p(x) over motif presence; [2].
-    p_y_given_x = JOINT / p_x[:, None]  # The likelihood p(y | x); [2, 2].
-    evidence = float(p_y_given_x[:, 1] @ p_x)  # p(y) = sum_x p(y | x) p(x), summed over the cause.
-    posterior = float(p_y_given_x[1, 1] * p_x[1]) / evidence  # p(X = 1 | Y = 1) = 0.36 / 0.44.
-    show("evidence p(Y = 1)", 0.44, evidence, quiet)
-    show("posterior p(X = 1 | Y = 1)", 0.818, posterior, quiet)
-    show("prior p(X = 1)", 0.60, float(p_x[1]), quiet)
-    show("likelihood ratio 0.6 / 0.2", 3.0, float(p_y_given_x[1, 1] / p_y_given_x[0, 1]), quiet)
-    return {"evidence": evidence, "posterior": posterior}
-
-
-# %% 11. Independence
-def independence(quiet: bool = False) -> dict:
-    """Section 1.11: p(x, y) equals p(x) p(y) under independence, and the table fails that test."""
-    p_x = JOINT.sum(dim=1)  # Marginal over motif presence; [2].
-    p_y = JOINT.sum(dim=0)  # Marginal over binding activity; [2].
-    product = float(p_x[1] * p_y[1])  # 0.60 times 0.44 = 0.264 under a factorized model.
-    joint_entry = float(JOINT[1, 1])  # The table's own entry, 0.36.
-    show("p(X = 1) p(Y = 1)", 0.264, product, quiet)
-    show("p(X = 1, Y = 1)", 0.36, joint_entry, quiet)
-    show("association missed by factorizing", 0.096, joint_entry - product, quiet)
-    return {"product": product, "joint_entry": joint_entry, "gap": joint_entry - product}
-
-
-# %% 12. The report
-def main() -> None:
-    parser = argparse.ArgumentParser(description=__doc__.split("\n")[0])
-    parser.add_argument("--samples", type=int, default=1000, help="largest Monte Carlo sample size")
-    parser.add_argument("--steps", type=int, default=5000, help="repeats of the ten-sample estimate")
-    parser.add_argument("--seed", type=int, default=0, help="random seed for the Monte Carlo draws")
-    parser.add_argument("--out", type=str, default="", help="directory for the run's output")
-    parser.add_argument("--quiet", action="store_true", help="print only the final report")
-    args = parser.parse_args()
-
-    set_seed(args.seed)  # Seed first, so the same seed gives the same Monte Carlo ladder.
-    if not args.quiet:
-        print("Chapter 1, Sections 1.3 to 1.11, worked example by worked example\n")
-    two_token_distribution(args.quiet)
-    continuous_densities(args.quiet)
-    adenine = expected_adenine(args.quiet)
-    estimates = monte_carlo(args.samples, args.steps, args.quiet)
-    spread = worked_variance(args.quiet)
-    peptides = peptide_covariance(args.quiet)
-    noising = noising_covariance(args.quiet)
-    table = joint_marginal_conditional(args.quiet)
-    posterior = bayes(args.quiet)
-    association = independence(args.quiet)
-
-    report = {
-        "E[h(X)], expected A count": adenine["expectation"],
-        "Monte Carlo estimate at N": estimates["ladder"][-1][0],
-        "Monte Carlo estimate": estimates["ladder"][-1][1],
-        "Var(h(X))": estimates["variance"],
-        "ten-sample averages taken": estimates["rounds"],
-        "sd of those averages": estimates["measured"],
-        "Var(X) on 0, 1, 2": spread["variance"],
-        "rho_XY across three peptides": peptides["rho"],
-        "Cov(X_t)_11 from diag(4, 1)": float(noising["anisotropic"][0, 0]),
-        "p(Y = 1), high binding": float(table["p_y"][1]),
-        "p(X = 1 | Y = 1)": posterior["posterior"],
-        "p(X=1,Y=1) - p(X=1)p(Y=1)": association["gap"],
-    }
-    print_report("Probability, Sections 1.3 to 1.11", report)
-    save_run(args.out, config=vars(args), report=report)
-
-
-if __name__ == "__main__":  # Execute when launched directly.
-    main()  # Work every example of Sections 1.3 to 1.11 and print the report.