Skip to main content

Week 10 — Which method would you choose?

Open in Colab · Download notebook · Week 10 textbook

English · Русский

Draft for instructor review. You can upload the downloaded notebook to Colab if the first link fails. Use a Python CPU runtime.

Run Setup, then work on the challenge announced in class. Predict a result, adapt the nearby examples to fill the marked lines, and inspect the images and table. Answer the Moodle poll for your own run and be ready to explain one failure.

Allow 5–10 minutes for setup and 15–20 minutes per challenge. The first three challenges compare reconstruction, retrieval and anomaly detection. Two optional challenges explore maps and generation. Compare quality and runtime within each experiment and decide what each metric leaves out.

TODO(instructor): confirm classroom pacing, activity selection and dataset before release.

Show plotting helpers 
# @title Plotting helpers (run once) { display-mode: "form" }

def image_rows(rows, caption):
    fig, axes = plt.subplots(len(rows), len(rows[0][1]), squeeze=False,
                             figsize=(len(rows[0][1]) * 1.15, len(rows) * 1.28))
    for row_axes, (name, images, labels) in zip(axes, rows):
        for i, (ax, pixels) in enumerate(zip(row_axes, images)):
            ax.imshow(pixels.reshape(8, 8), cmap="gray", vmin=0, vmax=1, interpolation="nearest")
            ax.set_xticks([])
            ax.set_yticks([])
            if i == 0:
                ax.set_ylabel(fill(name, width=15), rotation=0, ha="right", va="center", fontsize=8)
            if labels is not None:
                ax.set_title(str(labels[i]), fontsize=9)
    fig.suptitle(fill(caption, width=88), fontsize=10)
    plt.tight_layout(rect=(0, 0, 1, 0.92))
    plt.show()
    plt.close(fig)

def plot_maps(views, labels):
    fig, axes = plt.subplots(1, 3, figsize=(12, 4), layout="constrained")
    for ax, (name, coords) in zip(axes, views.items()):
        points = ax.scatter(*coords.T, c=labels, cmap="tab10", vmin=-0.5, vmax=9.5, s=10)
        ax.set_title(name)
        ax.set(xlabel="Coordinate 1", ylabel="Coordinate 2")
    fig.colorbar(points, ax=axes, ticks=range(10), label="Digit label")
    fig.suptitle("Same 500 images in all three maps. Color shows digit labels.")
    plt.show()
    plt.close(fig)

Setup

Run the folded plotting cell above, then run both setup cells.

Enter a birthday month and day without the year, or invent a valid date. Keep it in your notebook. Keep the printed variant for your own run. Moodle asks only which method worked best in each core task.

We use 1,797 offline 8×8 handwritten digits. The textbook experiments use 28×28 MNIST. A row of X has D = 64 pixels, and a row of H has d code coordinates. x_hat is a reconstruction, and y is a digit label. Labels help split the data, define the anomaly scenario and check results. Models fit pixel values.

The fixed split is approximately 60% training, 20% validation and 20% test. Keep the assigned settings for the test poll. Use validation for extra experiments.

from importlib.metadata import PackageNotFoundError, version
import platform
import shlex
import subprocess
import sys

packages = ("numpy", "scipy", "scikit-learn", "matplotlib", "torch", "threadpoolctl")

def package_versions():
    versions = {}
    for package in packages:
        try:
            versions[package] = version(package)
        except PackageNotFoundError:
            versions[package] = None
    return versions

before = package_versions()
missing = [name for name, installed in before.items() if installed is None]
if missing:
    command = [sys.executable, "-m", "pip", "install", *missing]
    if "google.colab" not in sys.modules:
        raise RuntimeError(f"Install missing packages and rerun Setup: {shlex.join(command)}")
    try:
        subprocess.run(command, check=True)
    except (OSError, subprocess.CalledProcessError) as exc:
        raise RuntimeError("Package setup failed. Try a fresh Colab session and rerun Setup.") from exc
    after = package_versions()
    loaded_names = {"scikit-learn": "sklearn"}
    if any(before[name] != after[name] and loaded_names.get(name, name) in sys.modules
           for name in packages):
        raise RuntimeError("Restart the session, then rerun Setup to load the installed versions.")
print("Python", platform.python_version(), "| CPU")
print(package_versions())
from datetime import date
from inspect import signature
import math
import random
from textwrap import fill
from time import perf_counter

import matplotlib.pyplot as plt
import numpy as np
import torch
from scipy.spatial.distance import cdist
from sklearn.cluster import KMeans
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
from sklearn.model_selection import train_test_split
from threadpoolctl import threadpool_limits
from torch import nn
from torch.utils.data import DataLoader, TensorDataset

torch.set_num_threads(1)
threadpool_limits(limits=1)

month = 3  # @param {type:"integer"}
day = 7  # @param {type:"integer"}

def variant_for(month, day):
    date(2000, month, day)
    seed = 100 * month + day
    variant = random.Random(seed).choice(range(3))
    return seed, "ABC"[variant], (4, 8, 16)[variant]

personal_seed, variant, d = variant_for(month, day)
print(f"Variant {variant} — {d} latent coordinates")

digits = load_digits()
X = digits.data.astype(np.float32) / 16
y = digits.target
train, held_out = train_test_split(
    np.arange(len(X)), test_size=0.4, stratify=y, random_state=1010
)
val, test = train_test_split(
    held_out, test_size=0.5, stratify=y[held_out], random_state=1010
)
assert X.shape == (1797, 64) and 0 <= X.min() <= X.max() <= 1
assert np.array_equal(np.sort(np.r_[train, val, test]), np.arange(len(X)))
print("Rows in train / validation / test:", len(train), len(val), len(test))

def tsne(perplexity, seed, iterations=750):
    # scikit-learn renamed n_iter to max_iter in 1.5.
    keyword = "max_iter" if "max_iter" in signature(TSNE).parameters else "n_iter"
    return TSNE(n_components=2, perplexity=perplexity, init="pca",
                learning_rate="auto", random_state=seed, **{keyword: iterations})

small = X[train[:24]]
assert PCA(2).fit_transform(small).shape == (24, 2)
assert KMeans(2, n_init=10, random_state=0).fit(small).cluster_centers_.shape == (2, 64)
assert np.isfinite(tsne(5, 0, 300).fit_transform(small)).all()
probe = torch.tensor([2.0], requires_grad=True)
probe.square().sum().backward()
assert probe.grad.item() == 4
print("CPU setup checks passed.")

Supplied helpers

Each activity uses a fresh model and task seed. Fit times are one observation. Query and score times are medians of 20 measurements after warm-up, alternating method order. Compare seconds on your own CPU. The tie thresholds are classroom conventions.

incomplete = "I could not obtain a complete result"

def run_activity(activity, offset):
    seed = personal_seed + offset
    torch.manual_seed(seed)
    np.random.seed(seed)
    try:
        activity(seed)
    except (AssertionError, TypeError, ValueError, IndexError, RuntimeError) as exc:
        print(f"{incomplete}. Check the marked lines. {exc}")

def timed(operation):
    start = perf_counter()
    result = operation()
    return result, perf_counter() - start

def batch_times(operations):
    samples = {name: [] for name in operations}
    for operation in operations.values():
        operation()
    names = list(operations)
    for repeat in range(20):
        for name in (names if repeat % 2 == 0 else names[::-1]):
            _, seconds = timed(operations[name])
            samples[name].append(seconds)
    return {name: float(np.median(times)) for name, times in samples.items()}

def report(scores, metric, question=None, tolerance=0, maximize=False, timing=False):
    assert np.isfinite(list(scores.values())).all(), "A comparison contains a nonfinite result."
    print(f"\n{metric}")
    for name, value in scores.items():
        print(f"{name:34} {value:.6g}")
    best = (max if maximize else min)(scores.values())
    margin = best * tolerance if timing or not maximize else tolerance
    winners = [name for name, value in scores.items() if abs(value - best) <= margin]
    tie = "Tie / too close to distinguish" if timing else "Tie / practically the same"
    prefix = f"Moodle Q{question}" if question is not None else "Result"
    print(f"{prefix}: {winners[0] if len(winners) == 1 else tie}")

class Autoencoder(nn.Module):
    def __init__(self, width):
        super().__init__()
        self.encoder = nn.Sequential(nn.Linear(64, 64), nn.ReLU(), nn.Linear(64, width))
        self.decoder = nn.Sequential(nn.Linear(width, 64), nn.ReLU(), nn.Linear(64, 64), nn.Sigmoid())

    def forward(self, x):
        return self.decoder(self.encoder(x))

def fit_network(model, dataset, seed, objective, epochs=40):
    batches = DataLoader(dataset, batch_size=128, shuffle=True,
                         generator=torch.Generator().manual_seed(seed))
    optimizer = torch.optim.Adam(model.parameters(), lr=0.001)
    for epoch in range(epochs):
        for batch_number, (x, target) in enumerate(batches):
            assert x.ndim == 2 and x.shape == target.shape and x.shape[1] == 64
            optimizer.zero_grad()
            loss = objective(model, x, target)
            assert torch.isfinite(loss), "Training loss is nonfinite."
            loss.backward()
            if epoch == 0 and batch_number == 0:
                gradients = [p.grad for p in model.parameters() if p.grad is not None]
                assert all(torch.isfinite(g).all() for g in gradients)
                assert any(torch.count_nonzero(g) for g in gradients), "No learning signal."
            optimizer.step()
    return model.eval()

def mse(model, x, target):
    return (model(x) - target).square().mean()

def clean_ae(images, seed):
    x = torch.from_numpy(images)
    return fit_network(Autoencoder(d), TensorDataset(x, x), seed, mse)

@torch.inference_mode()
def predict(operation, images):
    return operation(torch.from_numpy(images)).numpy()

def fit_pca(images, width):
    return PCA(n_components=width, svd_solver="full").fit(images)

def reconstruct_pca(model, images):
    return model.inverse_transform(model.transform(images))

Task 1 — Rescue a noisy digit

Predict: Which method will recover the strokes better? Which fits faster?

Complete the first convolution and the linear projection into the code. The shape path is [B, 1, 8, 8] → [B, 8, 8, 8] → [B, 512] → [B, d]. Use a 3×3 kernel with stride 1. The supplied dataset pairs noisy inputs with clean targets, and the decoder restores a flat 64-pixel image.

Compare against PCA fitted on clean training images. Both methods reconstruct the same clipped noisy test images. Score PCA’s raw output and use grayscale limits only for display. The noisy row is a reference.

Moodle Q1: Choose the lower clean-target MSE. Discuss fit time aloud. Errors within 1% and fit times within 10% of the minimum are treated as ties. Point to a digit where the numeric result and your visual judgment differ.

def denoising(seed):
    rng = np.random.default_rng(seed)
    noisy = np.clip(X + rng.normal(0, 0.35, X.shape), 0, 1).astype(np.float32)
    noisy_train, clean_train = map(torch.from_numpy, (noisy[train], X[train]))

    # Example: a 3×3 convolution from 2 to 4 channels preserves height and width.
    # nn.Conv2d(2, 4, kernel_size=3, stride=1, padding=1)
    convolution = None  # TODO(student): 1 input channel, 8 output channels
    # Example: nn.Linear(4 * 8 * 8, 6) maps a flattened 4-channel image to 6 coordinates.
    projection = None  # TODO(student): project the flattened features to d
    if convolution is None or projection is None:
        print(incomplete, "Complete the convolution and projection.")
        return
    with torch.no_grad():
        features = convolution(noisy_train[:2].reshape(2, 1, 8, 8))
        assert features.shape == (2, 8, 8, 8), "Check channels, stride and padding."
        assert projection(features.flatten(1)).shape == (2, d), "Check the linear dimensions."

    def fit_denoiser():
        model = Autoencoder(d)
        model.encoder = nn.Sequential(
            nn.Unflatten(1, (1, 8, 8)), convolution, nn.ReLU(), nn.Flatten(1), projection
        )
        return fit_network(model, TensorDataset(noisy_train, clean_train), seed, mse)

    pca, pca_seconds = timed(lambda: fit_pca(X[train], d))
    ae, ae_seconds = timed(fit_denoiser)
    assert predict(ae.encoder, X[test[:8]]).shape == (8, d)
    results = {
        "PCA": reconstruct_pca(pca, noisy[test]),
        "Denoising autoencoder": predict(ae, noisy[test]),
    }
    clean = X[test]
    print(f"Noisy image MSE (reference): {np.mean((noisy[test] - clean)**2):.6f}")
    report({name: np.mean((out - clean)**2) for name, out in results.items()},
           "Test MSE against clean images", 1, tolerance=0.01)
    report(dict(zip(results, (pca_seconds, ae_seconds))), "Fit seconds",
           tolerance=0.10, timing=True)
    rows = [("Clean", clean[:8], None), ("Noisy", noisy[test[:8]], None)]
    rows += [(name, out[:8], None) for name, out in results.items()]
    image_rows(rows, "Eight fixed test digits: clean target, noisy input and both reconstructions.")

run_activity(denoising, 100)

Task 2 — Find a matching digit

Predict: Does a useful reconstruction code also make a good search key?

Find one nearest training image for each test query. Compare original pixels, PCA codes and autoencoder codes using Euclidean distance in each space. Fit representations on clean training images. A match means the same digit label. The query and gallery have disjoint row IDs.

Select one gallery index per query, count equal digit labels, and turn the count into a top-1 accuracy. The distance matrix has one row per query and one column per gallery image.

Moodle Q2: Choose the most same-digit matches. Discuss query time aloud. A difference of one match or less counts as a tie. Timing includes encoding queries and searching the prepared gallery. Times within 10% of the minimum count as a tie.

def find_matches(encode, query, gallery):
    query_codes = encode(query)
    distances = cdist(query_codes, gallery, metric="sqeuclidean")
    assert distances.shape == (len(query), len(gallery))
    # Example: np.argmin(table, axis=0) selects a row for each column.
    nearest = None  # TODO(student): one minimum index per query
    if nearest is not None:
        assert nearest.shape == (len(query),) and ((0 <= nearest) & (nearest < len(gallery))).all()
    return nearest

def match_accuracy(found_labels, query_labels):
    # Example: int(np.sum(np.array([1, 4]) == np.array([1, 5]))) gives 1.
    count = None  # TODO(student): count equal labels
    if count is None:
        return None, None
    # Example: 3 correct predictions out of 5 gives 3 / 5.
    accuracy = None  # TODO(student): fraction of queries matched
    return count, accuracy

def retrieval(seed):
    probe = find_matches(np.asarray, X[test[:2]], X[train[:3]])
    count, accuracy = match_accuracy(np.array([1, 4]), np.array([1, 5]))
    if probe is None or count is None or accuracy is None:
        print(incomplete, "Complete the nearest-index and accuracy lines.")
        return
    assert (count, accuracy) == (1, 0.5), "Check label comparison and denominator."
    assert not np.intersect1d(train, test).size
    pca = fit_pca(X[train], d)
    ae = clean_ae(X[train], seed)
    encoders = {"Original pixels": np.asarray, "PCA codes": pca.transform,
                "Autoencoder codes": lambda x: predict(ae.encoder, x)}
    operations = {}
    for name, encode in encoders.items():
        gallery = encode(X[train])
        operations[name] = lambda encode=encode, gallery=gallery: find_matches(encode, X[test], gallery)
    matches = {name: train[operation()] for name, operation in operations.items()}
    scores = {}
    for name, ids in matches.items():
        count, accuracy = match_accuracy(y[ids], y[test])
        scores[name] = count
        print(f"{name}: {count}/{len(test)} ({accuracy:.1%})")
    report(scores, "Correct matches", 2, tolerance=1, maximize=True)
    report(batch_times(operations), "Query batch seconds", tolerance=0.10, timing=True)
    rows = [("Query", X[test[:8]], y[test[:8]])]
    rows += [(name, X[ids[:8]], y[ids[:8]]) for name, ids in matches.items()]
    image_rows(rows, "Eight test queries above their retrieved neighbors, with true labels.")

run_activity(retrieval, 200)

Task 3 — Find the intruder

Predict: Can an unfamiliar digit be easy to reconstruct?

Digits 0–8 define normal training data. In a mixed test batch of 80 normal images and 20 nines, flag ten images. Fit k-means, PCA and an autoencoder only on normal training rows. Labels define the task and reveal the result after scoring.

Complete the reduction from squared pixel residuals to one score per image. Then order the scores so the first rows are the most suspicious and count actual nines in the flagged set. At ranking time, original row IDs break ties.

Moodle Q3: Choose the most nines among the ten flags. Discuss scoring time aloud. Equal counts are a tie. Times within 10% of the minimum are a tie. The three methods’ raw scores are on different scales.

def reconstruction_scores(x, x_hat):
    squared = (x - x_hat)**2
    # Example: squared.mean(axis=0) gives one value per pixel across images.
    return None  # TODO(student): reduce the pixel axis

def rank_flags(scores):
    # Example: np.argsort(scores, kind="stable") orders the smallest scores first.
    return None  # TODO(student): descending stable argsort

def count_intruders(flag_labels):
    # Example: int(np.sum(flag_labels == 2)) counts the twos.
    return None  # TODO(student): count actual anomalies

def anomalies(seed):
    if reconstruction_scores(X[test[:2]], X[test[:2]]) is None or rank_flags(np.array([0., 1.])) is None or count_intruders(np.array([9, 2, 9])) is None:
        print(incomplete, "Complete the reduction, ranking and count.")
        return
    assert count_intruders(np.array([9, 2, 9])) == 2, "Count only the nines."
    rng = np.random.default_rng(seed)
    normal_train = train[y[train] != 9]
    normal_test, unfamiliar = test[y[test] != 9], test[y[test] == 9]
    ids = np.r_[rng.choice(normal_test, 80, replace=False),
                rng.choice(unfamiliar, 20, replace=False)]
    rng.shuffle(ids)
    assert len(np.unique(ids)) == 100 and (y[ids] == 9).sum() == 20
    assert (y[normal_train] != 9).all()
    # Scores use a shuffled batch. Original IDs break ties at the ranking step.
    x = X[ids]
    kmeans = KMeans(10, n_init=10, random_state=seed).fit(X[normal_train])
    pca = fit_pca(X[normal_train], d)
    ae = clean_ae(X[normal_train], seed)
    operations = {
        "Distance to k-means centers": lambda: cdist(x, kmeans.cluster_centers_, "sqeuclidean").min(axis=1),
        "PCA reconstruction error": lambda: reconstruction_scores(x, reconstruct_pca(pca, x)),
        "Autoencoder reconstruction error": lambda: reconstruction_scores(x, predict(ae, x)),
    }
    flags = {}
    row_order = np.argsort(ids)
    for name, operation in operations.items():
        scores = np.asarray(operation())
        assert scores.shape == (100,), "Average across pixels to get one score per image."
        scores = scores[row_order]
        rank = rank_flags(scores)
        assert np.array_equal(np.sort(rank), np.arange(100)), "Return 100 row indices."
        assert np.all(np.diff(scores[rank]) <= 0), "Put the highest scores first."
        flags[name] = ids[row_order[rank[:10]]]
    report({name: count_intruders(y[flag]) for name, flag in flags.items()},
           "Nines among ten flags (80 normal / 20 unfamiliar)", 3, maximize=True)
    report(batch_times(operations), "Score batch seconds, excluding ranking",
           tolerance=0.10, timing=True)
    image_rows([(name, X[flag], y[flag]) for name, flag in flags.items()],
               "The ten highest-scoring images per method with true labels revealed.")

run_activity(anomalies, 300)

Extension A — Can you trust this map?

Predict: Which map keeps more of each image’s ten closest neighbors?

Fit two-dimensional PCA and t-SNE maps to the same 500 images. Complete the step that excludes each image from its own neighbor list, then calculate neighbor retention. Use the diagonal of the square distance matrix. Labels color the maps after fitting.

Discuss: Compare PCA with perplexity-30 t-SNE. Report average original-neighbor retention and map construction time. Differences of at most 0.01 retention and times within 10% of the minimum are ties. Compare the extra perplexity-10 view visually. What could you infer about neighboring images? What would you avoid inferring about gaps between islands?

def map_neighbors(coordinates):
    distances = cdist(coordinates, coordinates, "sqeuclidean")
    # Example: setting a candidate distance to np.inf keeps it out of the nearest ten.
    self_distance = None  # TODO(student): value to put on the diagonal
    if self_distance is None:
        return None
    np.fill_diagonal(distances, self_distance)
    neighbors = np.argsort(distances, axis=1, kind="stable")[:, :10]
    assert not np.any(neighbors == np.arange(len(coordinates))[:, None]), "Exclude each row itself."
    return neighbors

def neighbor_overlap(original, mapped):
    # Example: np.intersect1d([0, 1], [1, 2]) returns array([1]).
    return None  # TODO(student): fraction of original neighbors retained

def maps(seed):
    ids, _ = train_test_split(train, train_size=500, stratify=y[train], random_state=1010)
    pixels = X[ids]
    original = map_neighbors(pixels)
    probe = neighbor_overlap(np.array([0, 1]), np.array([1, 2]))
    if original is None or probe is None:
        print(incomplete, "Complete self-exclusion and neighbor overlap.")
        return
    assert probe == 0.5, "One of two original neighbors is retained."
    pca_map, pca_seconds = timed(lambda: PCA(2, svd_solver="full").fit_transform(pixels))
    tsne_map, tsne_seconds = timed(lambda: tsne(30, seed).fit_transform(pixels))
    views = {"PCA map": pca_map, "t-SNE map (perplexity 30)": tsne_map}
    neighbors = {name: map_neighbors(coords) for name, coords in views.items()}
    retention = {
        name: np.mean([neighbor_overlap(a, b) for a, b in zip(original, mapped)])
        for name, mapped in neighbors.items()
    }
    report(retention, "Fraction of ten original neighbors retained", tolerance=0.01, maximize=True)
    report(dict(zip(views, (pca_seconds, tsne_seconds))), "Map construction seconds",
           tolerance=0.10, timing=True)
    views["t-SNE map (perplexity 10)"] = tsne(10, seed).fit_transform(pixels)
    plot_maps(views, y[ids])
    for i in (0, 1, 2):
        rows = [("Query + pixels", pixels[np.r_[i, original[i, :4]]], None)]
        rows += [(name, pixels[np.r_[i, near[i, :4]]], None) for name, near in neighbors.items()]
        image_rows(rows, f"Image {ids[i]} at left, followed by four neighbors per space.")

run_activity(maps, 400)

Extension B — Reconstruct or generate?

Predict: Will the better reconstructor also decode random codes into more recognizable digits?

Complete the sampled code and the sum of KL contributions for each image. The encoder predicts a mean and log variance. Convert log variance to a standard deviation and scale the supplied standard-normal noise before shifting by the mean.

The VAE minimizes Gaussian reconstruction negative log likelihood (pixel standard deviation 0.25) plus KL to a standard normal. Sum each term over its features, then average over the batch. The autoencoder trains with pixel MSE. Compare the two using the same test MSE, decoding the VAE encoder mean for a repeatable diagnostic.

Discuss: Choose lower test MSE (within 1% is a tie), then judge which random-code grid contains more recognizable digits. Both grids use the same 16 codes. Images show decoder means. The AE did not learn a standard-normal prior.

def sample_code(mu, logvar, epsilon):
    # Example: mean + std * noise shifts and scales standard-normal noise.
    std = None  # TODO(student): convert log variance to standard deviation
    return None if std is None else mu + std * epsilon

class VariationalAutoencoder(nn.Module):
    def __init__(self, width):
        super().__init__()
        self.encoder = nn.Sequential(nn.Linear(64, 64), nn.ReLU())
        self.mu = nn.Linear(64, width)
        self.logvar = nn.Linear(64, width)
        self.decoder = nn.Sequential(nn.Linear(width, 64), nn.ReLU(), nn.Linear(64, 64), nn.Sigmoid())

    def encode(self, x):
        h = self.encoder(x)
        return self.mu(h), self.logvar(h)

def kl_to_normal(mu, logvar):
    per_coordinate = 0.5 * (mu.square() + logvar.exp() - 1 - logvar)
    # Example: per_coordinate.sum(dim=0) gives one total per latent coordinate.
    return None  # TODO(student): sum the KL coordinates for each image

def vae_objective(model, x, target):
    mu, logvar = model.encode(x)
    z_sample = sample_code(mu, logvar, torch.randn_like(mu))
    assert z_sample.shape == mu.shape
    x_hat = model.decoder(z_sample)
    sd = 0.25
    nll = (0.5 * ((target - x_hat) / sd).square() + math.log(sd * math.sqrt(2 * math.pi))).sum(dim=1)
    return (nll + kl_to_normal(mu, logvar)).mean()

def generation(seed):
    zero, one = torch.zeros(2, d), torch.ones(2, d)
    sample = sample_code(zero, zero, one)
    if sample is None or kl_to_normal(zero, zero) is None:
        print(incomplete, "Complete the sampling and KL reduction.")
        return
    assert torch.allclose(sample, one), "At zero mean and unit variance the sample equals the noise."
    assert torch.allclose(sample_code(one, one * math.log(4), one), 3 * one), "Check the standard deviation."
    assert torch.equal(kl_to_normal(zero, zero), torch.zeros(2))
    assert torch.allclose(kl_to_normal(one, zero), torch.full((2,), d / 2)), "Check the sum axis."
    x = torch.from_numpy(X[train])
    ae = clean_ae(X[train], seed)
    vae = fit_network(VariationalAutoencoder(d), TensorDataset(x, x), seed, vae_objective)
    clean = X[test]
    outputs = {
        "Autoencoder": predict(ae, clean),
        "Variational autoencoder": predict(lambda x: vae.decoder(vae.encode(x)[0]), clean),
    }
    report({name: np.mean((out - clean)**2) for name, out in outputs.items()}, "Test MSE", tolerance=0.01)
    rows = [("Clean", clean[:8], None)] + [(name, out[:8], None) for name, out in outputs.items()]
    image_rows(rows, "Eight fixed test digits and their deterministic reconstructions.")
    codes = np.random.default_rng(seed).standard_normal((16, d)).astype(np.float32)
    grids = [("Autoencoder grid", predict(ae.decoder, codes)),
             ("Variational autoencoder grid", predict(vae.decoder, codes))]
    image_rows([(name if start == 0 else "", out[start:start + 8], None)
                for name, out in grids for start in (0, 8)],
               "Judge the two 16-image grids decoded from the same normal codes.")

run_activity(generation, 500)

Discuss your results

Choose a method for one job. Give one result that supports your choice, one image where it fails, and one cost you would consider before using it. Compare votes as observations from these runs, including ties and incomplete results.

TODO(instructor): confirm Moodle access and reporting, gradebook treatment, classroom timing and the latest Colab runtime before release.