Hands-on lab — Day 145: Overfitting and Underfitting

"""Overfitting and underfitting, decomposed and measured.

The two failures are usually described as a picture: a wiggly line through
every point, a straight line through none. That picture is true and it does
not tell you which one you have, because on your own data you cannot see
the true function.

What you can see is a decomposition. Every squared error splits into three
parts -- bias, variance and irreducible noise -- and the two failures are
one part each:

* **underfitting is bias**: the model class cannot represent the truth, so
  it is wrong in the same direction no matter what data you give it;
* **overfitting is variance**: the model class can represent the truth and
  much else besides, so it chases the noise and lands somewhere different
  for every training set;
* **noise is neither**, and it is the floor nothing can go below.

This module measures all three directly by fitting many models to many
independent training sets, and checks the decomposition against the error
that was actually observed. Everything is deterministic given a seed, and
everything runs on the CPU in seconds.
"""

from __future__ import annotations

import numpy as np

from sklearn.linear_model import LinearRegression, Ridge, SGDRegressor
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler

#: The function the data is actually generated from. A cubic, so a
#: degree-3 model can represent it exactly and a straight line cannot.
NOISE_SD = 2.0


def true_function(x):
    """The relationship the data really has: a cubic, and nothing else."""
    return 0.5 * x**3 - 2.0 * x + 1.0


def irreducible_variance() -> float:
    """The floor. No model of any capacity can score below this."""
    return NOISE_SD**2


def make_data(n: int, seed: int, noise_sd: float = NOISE_SD):
    """n points from the true cubic, plus independent Gaussian noise."""
    rng = np.random.default_rng(seed)
    x = np.sort(rng.uniform(-3.0, 3.0, n))
    y = true_function(x) + rng.normal(0.0, noise_sd, n)
    return x.reshape(-1, 1), y


def polynomial_model(degree: int, alpha: float = 0.0):
    """A polynomial fit of the given degree, optionally ridge-regularised.

    The `StandardScaler` between the features and the estimator is not
    decoration. Raw polynomial features up to degree 24 span many orders
    of magnitude and the normal equations become numerically hopeless;
    without scaling the *training* error starts rising with degree, which
    is a conditioning artefact masquerading as a result.
    """
    estimator = LinearRegression() if alpha == 0.0 else Ridge(alpha=alpha)
    return make_pipeline(PolynomialFeatures(degree), StandardScaler(), estimator)


def mse(model, X, y) -> float:
    """Mean squared error of a fitted model on the given rows."""
    return float(np.mean((model.predict(X) - y) ** 2))


# --------------------------------------------------------------------------
# 1. The capacity sweep
# --------------------------------------------------------------------------


def capacity_sweep(degrees, n_train: int = 25, train_seed: int = 145, test_seed: int = 246):
    """Train and test error at each polynomial degree.

    Rows are ``(degree, train_mse, test_mse, gap)``. Training error falls
    with capacity; test error does not, and where it turns is the whole
    subject.
    """
    X_train, y_train = make_data(n_train, train_seed)
    X_test, y_test = make_data(2000, test_seed)
    rows = []
    for degree in degrees:
        model = polynomial_model(degree).fit(X_train, y_train)
        train = mse(model, X_train, y_train)
        test = mse(model, X_test, y_test)
        rows.append((degree, round(train, 4), round(test, 4), round(test - train, 4)))
    return rows


def best_degree(rows) -> int:
    """The degree with the lowest test error in a capacity sweep."""
    return min(rows, key=lambda row: row[2])[0]


# --------------------------------------------------------------------------
# 2. Regularisation
# --------------------------------------------------------------------------


def regularisation_sweep(alphas, degree: int = 24, n_train: int = 25):
    """Train and test error at each ridge penalty, holding capacity fixed.

    The model class does not change here. Only how much the fit is
    discouraged from using it does -- which is the point: overfitting is
    not a property of the model class alone.
    """
    X_train, y_train = make_data(n_train, 145)
    X_test, y_test = make_data(2000, 246)
    rows = []
    for alpha in alphas:
        model = polynomial_model(degree, alpha=alpha).fit(X_train, y_train)
        rows.append((alpha, round(mse(model, X_train, y_train), 4), round(mse(model, X_test, y_test), 4)))
    return rows


def best_alpha(rows) -> float:
    """The penalty with the lowest test error."""
    return min(rows, key=lambda row: row[2])[0]


# --------------------------------------------------------------------------
# 3. What more data fixes, and what it does not
# --------------------------------------------------------------------------


def data_sweep(sizes, degrees=(1, 4, 24), test_seed: int = 246):
    """Test error against training-set size, at three capacities.

    Rows are ``(n, {degree: test_mse})``. The three columns behave
    completely differently, and the difference is the practical content of
    the whole lesson.
    """
    X_test, y_test = make_data(2000, test_seed)
    rows = []
    for n in sizes:
        X_train, y_train = make_data(n, 900 + n)
        scores = {}
        for degree in degrees:
            model = polynomial_model(degree).fit(X_train, y_train)
            scores[degree] = round(mse(model, X_test, y_test), 4)
        rows.append((n, scores))
    return rows


# --------------------------------------------------------------------------
# 4. The decomposition itself
# --------------------------------------------------------------------------


def bias_variance(degree: int, n_train: int = 25, datasets: int = 200, grid: int = 200):
    """Measure bias squared, variance and noise for one model class.

    Fit ``datasets`` models, each to its own independent training set, and
    predict the same fixed grid of query points with all of them. Then:

    * **bias squared** is how far the *average* prediction sits from the
      truth -- an error the model class makes every time;
    * **variance** is how much the predictions scatter around their own
      average -- an error that changes with the training set;
    * **noise** is the irreducible term, known here because we generated it.

    Returns a dict with all three, their sum, and the squared error
    actually observed against freshly noisy targets. Those last two should
    agree, and checking that they do is what makes this a measurement
    rather than a recitation.
    """
    query = np.linspace(-3.0, 3.0, grid).reshape(-1, 1)
    truth = true_function(query.ravel())

    predictions = np.empty((datasets, grid))
    for i in range(datasets):
        X_train, y_train = make_data(n_train, 10_000 + i)
        model = polynomial_model(degree).fit(X_train, y_train)
        predictions[i] = model.predict(query)

    mean_prediction = predictions.mean(axis=0)
    bias_squared = float(np.mean((mean_prediction - truth) ** 2))
    variance = float(np.mean(predictions.var(axis=0)))
    noise = irreducible_variance()

    observed_rng = np.random.default_rng(7)
    noisy_targets = truth + observed_rng.normal(0.0, NOISE_SD, predictions.shape)
    observed = float(np.mean((predictions - noisy_targets) ** 2))

    return {
        "bias_squared": round(bias_squared, 4),
        "variance": round(variance, 4),
        "noise": round(noise, 4),
        "predicted_total": round(bias_squared + variance + noise, 4),
        "observed": round(observed, 4),
    }


def decomposition_table(degrees, n_train: int = 25, datasets: int = 200):
    """The decomposition at several capacities, as one table."""
    return [(degree, bias_variance(degree, n_train, datasets)) for degree in degrees]


# --------------------------------------------------------------------------
# 5. Early stopping
# --------------------------------------------------------------------------


def training_history(epochs: int = 600, degree: int = 14, n_train: int = 25, seed: int = 0):
    """Train and test error after every epoch of gradient descent.

    Capacity is fixed and the data is fixed. The only thing that changes is
    how long the fit has been allowed to run -- which turns out to be a
    capacity knob of its own.
    """
    X_train, y_train = make_data(n_train, 145)
    X_test, y_test = make_data(2000, 246)

    features = PolynomialFeatures(degree)
    scaler = StandardScaler()
    P_train = scaler.fit_transform(features.fit_transform(X_train))
    P_test = scaler.transform(features.transform(X_test))

    model = SGDRegressor(
        learning_rate="constant", eta0=0.003, penalty=None, random_state=seed
    )
    train_history, test_history = [], []
    for _epoch in range(epochs):
        model.partial_fit(P_train, y_train)
        train_history.append(float(np.mean((model.predict(P_train) - y_train) ** 2)))
        test_history.append(float(np.mean((model.predict(P_test) - y_test) ** 2)))
    return train_history, test_history


def is_monotonically_decreasing(values, tolerance: float = 1e-9) -> bool:
    """Whether a sequence never increases, within a numerical tolerance."""
    return all(a >= b - tolerance for a, b in zip(values, values[1:]))


def first_increase(values) -> int:
    """The index of the first value larger than the one before it."""
    for i, (a, b) in enumerate(zip(values, values[1:]), start=1):
        if b > a:
            return i
    return len(values)


def stop_with_patience(values, patience: int) -> int:
    """Index of the epoch a patience-based early stop would have chosen.

    Stops once ``patience`` consecutive epochs have failed to improve on
    the best seen so far, and returns the best epoch rather than the epoch
    it stopped at -- which is what a real implementation restores.
    """
    best_index, best_value, waited = 0, values[0], 0
    for i, value in enumerate(values[1:], start=1):
        if value < best_value:
            best_index, best_value, waited = i, value, 0
        else:
            waited += 1
            if waited >= patience:
                break
    return best_index

#!/usr/bin/env python3
"""Print every measured pair in this lab as one table.

The harness compares this output byte for byte against
expected-output/measured-values.txt, so the report is not a convenience:
it is how the lab notices that a number in the lesson has gone stale.
"""

import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).resolve().parent))

import numpy as np  # noqa: E402

import fitting_lib as f  # noqa: E402


def rule(title: str) -> None:
    print()
    print(title)
    print("-" * len(title))


def main() -> None:
    print("Day 145 -- overfitting and underfitting, decomposed and measured")
    print("=" * 63)
    print()
    print(f"  The data is a cubic plus Gaussian noise of standard deviation {f.NOISE_SD}.")
    print(f"  So no model of any capacity can score below {f.irreducible_variance():.4f}.")

    rule("1. The capacity sweep: 25 training rows, 2000 test rows")
    degrees = [1, 2, 3, 4, 6, 8, 10, 14, 18, 24]
    rows = f.capacity_sweep(degrees)
    print("   deg    train MSE        test MSE             gap")
    for degree, train, test, gap in rows:
        print(f"  {degree:4d}   {train:10.4f}   {test:13.4f}   {gap:13.4f}")
    print(f"  lowest test error at degree {f.best_degree(rows)}")
    train_column = [row[1] for row in rows]
    print(f"  training error falls at every step through degree 14 : "
          f"{f.is_monotonically_decreasing(train_column[:8])}")
    print("  after that it wobbles by 0.0636, which is numerical rather than statistical:")
    print("  25 training rows and degree 24 supplies exactly 25 features")

    rule("2. The same degree-24 model, with a penalty")
    reg = f.regularisation_sweep([0.0, 1e-6, 1e-3, 0.1, 1.0, 10.0, 100.0])
    print("      alpha    train MSE         test MSE")
    for alpha, train, test in reg:
        print(f"  {alpha:>9}   {train:10.4f}   {test:14.4f}")
    print(f"  lowest test error at alpha {f.best_alpha(reg)}")
    print(f"  the penalty improves the test error by a factor of "
          f"{reg[0][2] / min(r[2] for r in reg):.0f}")
    print("  training error rises with every increase in the penalty; that is the trade")

    rule("3. What more data fixes, and what it does not")
    data = f.data_sweep([15, 25, 50, 100, 400, 2000])
    print("      n    degree 1     degree 4        degree 24")
    for n, scores in data:
        print(f"  {n:5d}   {scores[1]:9.4f}   {scores[4]:9.4f}   {scores[24]:14.4f}")
    underfit = [scores[1] for _n, scores in data]
    overfit = [scores[24] for _n, scores in data]
    print(f"  underfit model, 15 rows to 2000 : {underfit[0]:.4f} -> {underfit[-1]:.4f}")
    print(f"  overfit model,  15 rows to 2000 : {overfit[0]:.4f} -> {overfit[-1]:.4f}")
    print(f"  the irreducible floor           : {f.irreducible_variance():.4f}")
    print("  more data cures overfitting completely and underfitting not at all")
    print(f"  the degree-24 column peaks at n=25, where features ({25}) equal rows")

    rule("4. The decomposition: bias squared, variance, noise")
    print("   deg     bias^2       variance      noise      predicted     observed")
    for degree, result in f.decomposition_table([1, 2, 3, 4, 6, 8, 12]):
        print(
            f"  {degree:4d}  {result['bias_squared']:10.4f}  {result['variance']:13.4f}  "
            f"{result['noise']:9.4f}  {result['predicted_total']:13.4f}  {result['observed']:13.4f}"
        )
    print("  underfitting is bias; overfitting is variance; the sum is the error")
    print("  degree 2 has MORE bias than degree 1 and twice the variance:")
    print("  the true function is odd, so a quadratic term buys nothing and costs")

    rule("5. Early stopping: the same model, trained for longer")
    train_history, test_history = f.training_history()
    best = int(np.argmin(test_history))
    print(f"  epochs run                        : {len(train_history)}")
    print(f"  training error falls every epoch  : {f.is_monotonically_decreasing(train_history)}")
    print(f"  training error, first to last     : {train_history[0]:.4f} -> {train_history[-1]:.4f}")
    print(f"  best test error                   : {test_history[best]:.4f} at epoch {best + 1}")
    print(f"  test error at epoch 600           : {test_history[-1]:.4f}")
    print(f"  worst test error after the best   : {max(test_history[best + 1:]):.4f}")
    print(f"  epochs worse than the best        : {sum(1 for v in test_history if v > min(test_history))} of 600")
    print(f"  generalisation gap, epoch 1       : {test_history[0] - train_history[0]:.4f}")
    print(f"  generalisation gap, epoch 600     : {test_history[-1] - train_history[-1]:.4f}")
    for patience in (5, 10, 20, 50):
        chosen = f.stop_with_patience(test_history, patience)
        print(f"  patience {patience:<3} would stop at epoch  : {chosen + 1} (test {test_history[chosen]:.4f})")
    print("  the test curve is not a clean U: it rises, then partly recovers,")
    print("  without ever again beating the value it reached at epoch 14")


if __name__ == "__main__":
    main()

"""The reference solutions: the two failures, decomposed and measured.

Every number here was captured from a real run of this file on the
authoring machine. If a number changes, the claim in the lesson is wrong
and one of the two must be fixed.
"""

import numpy as np
import pytest

import fitting_lib as f


# --- 1. The capacity sweep -----------------------------------------------


def test_01_training_error_falls_with_capacity_and_test_error_does_not():
    rows = f.capacity_sweep([1, 2, 3, 4, 6, 8, 10, 14, 18, 24])
    assert rows == [
        (1, 11.3217, 9.8274, -1.4942),
        (2, 11.3173, 9.8801, -1.4372),
        (3, 2.7076, 6.123, 3.4154),
        (4, 2.4964, 5.4911, 2.9948),
        (6, 1.9569, 15.8217, 13.8648),
        (8, 1.701, 26.1708, 24.4697),
        (10, 1.357, 528.4798, 527.1227),
        (14, 0.9685, 31307.2782, 31306.3097),
        (18, 1.0037, 75539.3618, 75538.3581),
        (24, 1.0321, 226667.4689, 226666.4368),
    ]
    # Test error is U-shaped: it falls, bottoms, then explodes.
    assert f.best_degree(rows) == 4
    test = [row[2] for row in rows]
    assert test[0] > test[3] < test[-1]
    assert test[-1] / test[3] > 40_000


def test_01b_training_error_falls_until_the_numerics_give_out():
    rows = f.capacity_sweep([1, 2, 3, 4, 6, 8, 10, 14, 18, 24])
    train = [row[1] for row in rows]
    # Monotone through degree 14, then it wobbles by about 0.06.
    assert f.is_monotonically_decreasing(train[:8])
    assert not f.is_monotonically_decreasing(train)
    assert round(max(train[7:]) - min(train[7:]), 4) == 0.0636
    # That wobble is numerical, not statistical: the training set has 25
    # rows and degree 24 supplies exactly 25 features.
    assert train[7] < 1.0 < train[-1]


def test_01c_the_generalisation_gap_is_the_diagnostic():
    rows = f.capacity_sweep([1, 2, 3, 4, 6, 8, 10, 14, 18, 24])
    gaps = {degree: gap for degree, _tr, _te, gap in rows}
    # Underfitting: the gap is NEGATIVE. Test error is below training
    # error, because the model is too rigid to have chased any noise.
    assert gaps[1] < 0 and gaps[2] < 0
    # The best model has a small positive gap.
    assert 0 < gaps[4] < 3.0
    # Overfitting: the gap is the error.
    assert gaps[24] > 200_000


# --- 2. Regularisation ---------------------------------------------------


def test_02_a_penalty_rescues_the_same_model_class():
    rows = f.regularisation_sweep([0.0, 1e-6, 1e-3, 0.1, 1.0, 10.0, 100.0])
    assert rows == [
        (0.0, 1.0321, 226667.4689),
        (1e-06, 1.2031, 130776.6548),
        (0.001, 1.5741, 128.3127),
        (0.1, 2.2259, 15.8339),
        (1.0, 2.7461, 5.7257),
        (10.0, 3.9689, 6.1559),
        (100.0, 6.28, 6.784),
    ]
    # Same degree, same data. Only the penalty changed.
    assert f.best_alpha(rows) == 1.0
    unpenalised = rows[0][2]
    best = min(row[2] for row in rows)
    assert round(unpenalised / best, 0) == 39588.0


def test_02b_the_penalty_trades_training_error_for_test_error():
    rows = f.regularisation_sweep([0.0, 1e-6, 1e-3, 0.1, 1.0, 10.0, 100.0])
    train = [row[1] for row in rows]
    # Training error rises monotonically with the penalty, always.
    assert all(a < b for a, b in zip(train, train[1:]))
    # Test error is U-shaped in the penalty too: too much is also wrong.
    test = [row[2] for row in rows]
    assert test[4] < test[5] < test[6]
    # And the best-regularised degree-24 model is close to the best
    # unregularised degree-4 one, from a completely different direction.
    sweep = f.capacity_sweep([4])
    assert abs(min(test) - sweep[0][2]) < 0.3


# --- 3. What more data fixes ---------------------------------------------


def test_03_more_data_cures_overfitting_and_does_nothing_for_underfitting():
    rows = f.data_sweep([15, 25, 50, 100, 400, 2000])
    assert rows == [
        (15, {1: 8.5023, 4: 4.9218, 24: 215413.2388}),
        (25, {1: 8.862, 4: 6.1904, 24: 64631547.2994}),
        (50, {1: 8.2457, 4: 4.2661, 24: 6070.3302}),
        (100, {1: 8.3583, 4: 4.2934, 24: 5.3571}),
        (400, {1: 8.3007, 4: 3.9958, 24: 4.3139}),
        (2000, {1: 8.2393, 4: 3.988, 24: 4.0055}),
    ]
    underfit = [scores[1] for _n, scores in rows]
    overfit = [scores[24] for _n, scores in rows]
    # The underfit model is flat: 133 times more data buys 0.26.
    assert round(max(underfit) - min(underfit), 4) == 0.6227
    assert abs(underfit[0] - underfit[-1]) < 0.3
    # The overfit model falls by seven orders of magnitude.
    assert overfit[1] / overfit[-1] > 1e7


def test_03b_both_good_models_converge_to_the_irreducible_floor():
    rows = f.data_sweep([15, 25, 50, 100, 400, 2000])
    floor = f.irreducible_variance()
    assert floor == 4.0
    at_2000 = dict(rows[-1][1])
    # Degree 4 and degree 24 both land on the floor, from opposite sides.
    assert abs(at_2000[4] - floor) < 0.02
    assert abs(at_2000[24] - floor) < 0.01
    # The underfit model never gets near it, at any amount of data.
    assert at_2000[1] > floor * 2


def test_03c_the_overfit_column_peaks_at_the_interpolation_threshold():
    """Degree 24 supplies exactly 25 features, so n=25 is the worst case."""
    rows = f.data_sweep([15, 25, 50, 100, 400, 2000])
    overfit = [scores[24] for _n, scores in rows]
    # Not monotone: it gets worse before it gets better.
    assert overfit[1] > overfit[0]
    assert overfit[1] == max(overfit)
    # And the peak is exactly where features equal rows.
    from sklearn.preprocessing import PolynomialFeatures

    n_features = PolynomialFeatures(24).fit_transform(np.zeros((3, 1))).shape[1]
    assert n_features == 25
    assert rows[1][0] == n_features


# --- 4. The decomposition ------------------------------------------------


def test_04_underfitting_is_bias_and_overfitting_is_variance():
    underfit = f.bias_variance(1)
    right = f.bias_variance(3)
    overfit = f.bias_variance(12)

    assert underfit["bias_squared"] == 4.2985
    assert underfit["variance"] == 0.7112
    assert right["bias_squared"] == 0.0033
    assert right["variance"] == 0.8399
    assert overfit["bias_squared"] == 2803.5354
    assert overfit["variance"] == 452183.1336

    # Underfitting: bias dominates variance by a factor of six.
    assert underfit["bias_squared"] / underfit["variance"] > 6
    # The true function is a cubic, so at degree 3 the bias vanishes.
    assert right["bias_squared"] < 0.01
    # Overfitting: variance dominates bias by a factor of a hundred.
    assert overfit["variance"] / overfit["bias_squared"] > 100


def test_04b_the_decomposition_predicts_the_error_that_was_observed():
    """Bias squared plus variance plus noise equals the error, measured.

    The tolerance is 1.1 percent rather than something tighter because the
    observed term is itself a Monte Carlo estimate over 40,000 noisy
    targets and carries its own sampling error. The worst disagreement
    across seven capacities is 1.003 percent, at degree 6; five of the
    seven agree to better than a quarter of a percent.
    """
    worst = 0.0
    for degree in (1, 2, 3, 4, 6, 8, 12):
        result = f.bias_variance(degree)
        predicted = result["predicted_total"]
        observed = result["observed"]
        relative = abs(predicted - observed) / observed
        worst = max(worst, relative)
        assert relative < 0.011, (degree, result)
        # And the three parts sum to the total. The tolerance is 0.0002
        # because each part is stored already rounded to four places, so
        # summing the rounded parts can differ from the rounded sum by up
        # to half a unit in the last place per part.
        parts = result["bias_squared"] + result["variance"] + result["noise"]
        assert abs(parts - predicted) <= 0.0002, (degree, parts, predicted)
    assert round(worst, 5) == 0.01003


def test_04c_a_bigger_model_class_is_not_always_less_biased():
    """Degree 2 has MORE bias than degree 1, on an odd true function."""
    one = f.bias_variance(1)
    two = f.bias_variance(2)
    assert one["bias_squared"] == 4.2985
    assert two["bias_squared"] == 4.3342
    assert two["bias_squared"] > one["bias_squared"]
    # And it pays for that with double the variance.
    assert two["variance"] > 2 * one["variance"] - 0.01
    # So degree 2 is worse than degree 1 on both counts.
    assert two["predicted_total"] > one["predicted_total"]


# --- 5. Early stopping ---------------------------------------------------


def test_05_training_longer_makes_training_error_better_and_the_model_worse():
    train, test = f.training_history()
    assert len(train) == len(test) == 600
    # Training error falls at every single epoch, for six hundred epochs.
    assert f.is_monotonically_decreasing(train)
    assert round(train[0], 4) == 7.3906
    assert round(train[-1], 4) == 2.4744
    # Test error bottoms at epoch 14 and is worse at 600.
    best = int(np.argmin(test))
    assert best == 13
    assert round(test[best], 4) == 5.4555
    assert round(test[-1], 4) == 5.8978
    assert test[-1] > test[best]


def test_05b_the_generalisation_gap_grows_while_training_error_falls():
    train, test = f.training_history()
    first_gap = test[0] - train[0]
    last_gap = test[-1] - train[-1]
    assert round(first_gap, 4) == 0.6771
    assert round(last_gap, 4) == 3.4234
    assert last_gap / first_gap > 5


def test_05c_the_test_curve_is_not_a_clean_u_which_is_why_patience_exists():
    _train, test = f.training_history()
    best = int(np.argmin(test))
    after = test[best + 1 :]
    # It rises to 7.1435, then recovers to 5.8978 -- without ever again
    # beating the 5.4555 it reached at epoch 14.
    assert round(max(after), 4) == 7.1435
    assert round(test[-1], 4) == 5.8978
    assert min(after) > test[best]
    # 599 of the 600 epochs are worse than the best one.
    assert sum(1 for v in test if v > min(test)) == 599
    # Here every patience from 5 to 50 recovers the true best epoch --
    # but that is luck on this run, not a property of the rule.
    for patience in (5, 10, 20, 50):
        assert f.stop_with_patience(test, patience) == best
    assert f.first_increase(test) == best + 1

"""Machinery checks: the helpers behave, before any claim is made.

These four tests are solved in both `starter/` and `examples/`. They exist
so that a broken helper reports itself as a broken helper rather than as a
surprising scientific result.
"""

import numpy as np
import pytest

import fitting_lib as f


def test_the_data_generator_is_the_true_function_plus_noise():
    X, y = f.make_data(2000, seed=1)
    residuals = y - f.true_function(X.ravel())
    # The noise is centred and has the standard deviation it claims.
    assert abs(float(residuals.mean())) < 0.15
    assert abs(float(residuals.std()) - f.NOISE_SD) < 0.1
    assert f.irreducible_variance() == f.NOISE_SD**2 == 4.0
    # Two calls at one seed agree; two seeds do not.
    assert np.array_equal(f.make_data(50, 3)[1], f.make_data(50, 3)[1])
    assert not np.array_equal(f.make_data(50, 3)[1], f.make_data(50, 4)[1])


def test_a_degree_three_model_can_represent_the_truth_exactly():
    # With no noise at all, a cubic fit should be essentially perfect and
    # a straight line should not.
    X, y = f.make_data(200, seed=2, noise_sd=0.0)
    cubic = f.polynomial_model(3).fit(X, y)
    line = f.polynomial_model(1).fit(X, y)
    assert f.mse(cubic, X, y) < 1e-12
    # A straight line cannot, and is left with 4.174 of pure bias.
    assert round(f.mse(line, X, y), 3) == 4.174


def test_scaling_is_what_keeps_the_high_degree_fit_numerically_sane():
    from sklearn.linear_model import LinearRegression
    from sklearn.pipeline import make_pipeline
    from sklearn.preprocessing import PolynomialFeatures

    X, y = f.make_data(25, seed=145)
    scaled = f.polynomial_model(24).fit(X, y)
    unscaled = make_pipeline(PolynomialFeatures(24), LinearRegression()).fit(X, y)
    # Both fit the training data; the scaled pipeline fits it better,
    # because the unscaled normal equations are badly conditioned.
    assert f.mse(scaled, X, y) < f.mse(unscaled, X, y)


def test_the_stopping_helpers_agree_on_a_hand_checkable_sequence():
    values = [10.0, 8.0, 5.0, 6.0, 7.0, 9.0]
    assert f.first_increase(values) == 3
    assert not f.is_monotonically_decreasing(values)
    assert f.is_monotonically_decreasing([5.0, 4.0, 4.0, 3.0])
    # Patience 2 waits two non-improving epochs, then returns the best.
    assert f.stop_with_patience(values, 2) == 2
    # A sequence that keeps improving is never stopped early.
    assert f.stop_with_patience([5.0, 4.0, 3.0, 2.0], 2) == 3

lesson_id: D145
day: 145
kind: guided-build
languages:
  - python
  - bash
setup_commands:
  - cd labs/sections/machine-learning/day-145-overfitting-and-underfitting
  - python3 -m venv .venv
  - .venv/bin/pip install -r requirements/requirements.txt
  - >-
    .venv/bin/python3 -c "import numpy, sklearn; print(numpy.__version__,
    sklearn.__version__)"
run_commands:
  - .venv/bin/pytest examples -q
  - .venv/bin/pytest starter -q
  - .venv/bin/python3 examples/report_measurements.py
test_commands:
  - bash tests/run_tests.sh
cleanup_commands:
  - >-
    find . -path ./.venv -prune -o -type d -name '__pycache__' -print -exec rm -rf -- {}
    +
  - rm -rf .pytest_cache
  - 'rm -rf .venv  # optional: removes the lab virtual environment'
  - 'git checkout -- starter/  # optional: reset your work'
requires_network: true
requires_api_key: false
estimated_minutes: 60
last_executed: '2026-08-27'
executed_on: >-
  macOS 26.5.2 (Apple Silicon, arm64, CPU only -- no GPU is needed or used), Python
  3.14.0, numpy 2.5.2, scikit-learn 1.9.0, pytest 9.1.1, bash 3.2.57 -- bash
  tests/run_tests.sh -> 14 checks, 0 failure(s), exit 0. pytest examples -q -> 18 passed.
  pytest starter -q -> 4 passed, 14 skipped (the four machinery checks in
  test_fitting_lib.py are solved in both directories; the fourteen exercise stubs in
  starter/test_fitting_claims.py are untouched). Everything ran through a real lab-local
  .venv created by the documented setup commands; scikit-learn pulled in scipy 1.18.1,
  joblib 1.5.3 and threadpoolctl 3.6.0 as its own dependencies, none of which this lab
  imports directly. The lab is fully offline after the pip install -- every dataset is
  generated on the spot from a seeded numpy.random.default_rng, nothing is downloaded, no
  dataset is bundled, and harness check 9 confirms no URL appears anywhere in starter/ or
  examples/ source. Section 7 of the harness copies examples/ into a mktemp-d scratch
  directory, confirms 18 passed, rewrites `assert f.best_degree(rows) == 4` to 24,
  confirms a non-zero exit naming the failing test, and removes the scratch directory.
  Separately, by hand, `assert underfit["bias_squared"] == 4.2985` was changed to 0.0001
  in examples/test_fitting_claims.py and the whole harness re-run: it reported 14 checks,
  2 failure(s) and exited 1 (both the pytest run and the pytest-free direct reproduction
  in section 2 caught it); the file was restored and the harness returned to 14 checks, 0
  failure(s), exit 0. MEASURED PAIRS, all captured verbatim in
  expected-output/measured-values.txt. The data is a cubic plus Gaussian noise of standard
  deviation 2.0, so the irreducible floor is exactly 4.0000. (1) THE CAPACITY SWEEP on 25
  training rows and 2000 test rows: training MSE falls 11.3217, 11.3173, 2.7076, 2.4964,
  1.9569, 1.7010, 1.3570, 0.9685 at degrees 1, 2, 3, 4, 6, 8, 10, 14 while test MSE runs
  9.8274, 9.8801, 6.1230, 5.4911, 15.8217, 26.1708, 528.4798, 31307.2782 and on to
  226667.4689 at degree 24 -- a U-curve whose right-hand arm is more than 40000 times its
  minimum. (2) THE SIGN OF THE GAP IS THE DIAGNOSTIC: at degrees 1 and 2 the gap is
  NEGATIVE, -1.4942 and -1.4372, because a model too rigid to chase noise has none to be
  flattered by; this is the signature of underfitting and the one most often mistaken for
  a broken split. (3) REGULARISATION: the same degree-24 model with ridge penalties 0.0,
  1e-6, 1e-3, 0.1, 1.0, 10.0 and 100.0 scores test MSE 226667.4689, 130776.6548, 128.3127,
  15.8339, 5.7257, 6.1559 and 6.7840 while its training MSE rises monotonically 1.0321,
  1.2031, 1.5741, 2.2259, 2.7461, 3.9689, 6.2800 -- an improvement of a factor of 39588,
  bought precisely by fitting the training data less well, and with its own U-curve in the
  penalty. The best-regularised degree-24 model lands within 0.3 of the best unregularised
  degree-4 model, the same destination from a different direction. (4) MORE DATA: from 15
  rows to 2000 the degree-1 model moves 8.5023 -> 8.2393, a total range of 0.6227 across a
  133-fold increase, while the degree-24 model moves 215413.2388 -> 4.0055, seven orders of
  magnitude, landing exactly on the irreducible floor; degree 4 lands at 3.9880. More data
  cures variance completely and bias not at all. (5) THE DECOMPOSITION, measured over 200
  independent training sets: at degree 1 bias squared 4.2985 against variance 0.7112; at
  degree 3 bias 0.0033 against variance 0.8399 (the true function IS a cubic); at degree 12
  bias 2803.5354 against variance 452183.1336. The predicted totals 9.0097, 4.8432 and
  454990.6691 match the separately observed squared errors 9.0295, 4.8431 and 455027.8625.
  (6) EARLY STOPPING: over 600 epochs of gradient descent on fixed data and a fixed model,
  training MSE falls at EVERY epoch from 7.3906 to 2.4744 while test MSE bottoms at 5.4555
  at epoch 14 and is 5.8978 at epoch 600; the generalisation gap grows from 0.6771 to
  3.4234. Training time is a capacity dial. FIVE HONESTY CALLS. FIRST: the training column
  is NOT monotone over the whole sweep -- it wobbles upward by 0.0636 past degree 14 -- and
  the lab asserts monotonicity through degree 14 while asserting that it fails overall,
  rather than pretending the curve is clean. The cause is numerical: 25 rows and exactly 25
  polynomial features make the system square and catastrophically ill-conditioned. SECOND,
  and found by accident: the degree-24 column of the data sweep is NOT monotone in n. It is
  worse at 25 rows (64631547.2994) than at 15 (215413.2388), because 25 is exactly the
  interpolation threshold where features equal rows; at 15 rows least squares returns the
  minimum-norm solution, which is quietly a form of regularisation. The lab proves this
  with a feature count rather than asserting it, and the lesson names it as the left-hand
  edge of double descent. THIRD: degree 2 measures MORE bias than degree 1 (4.3342 against
  4.2985) despite containing it as a special case, plus double the variance. The true
  function is odd so a quadratic term buys nothing and still costs; the lesson states the
  mechanism rather than claiming a theorem, and notes that in population terms a larger
  class cannot have more bias. FOURTH: the decomposition's agreement with the observed
  error is asserted at 1.1 percent rather than something tighter, because the observed
  column is itself a Monte Carlo estimate; the worst disagreement is 0.01003, at degree 6,
  and five of seven capacities agree to better than 0.25 percent. Separately the three
  parts are asserted to sum to the total within 0.0002 rather than exactly, because each
  part is stored already rounded. FIFTH: the test curve after its minimum is not a clean U
  -- it rises to 7.1435 around epoch 84 and partly recovers to 5.8978 -- and although every
  patience from 5 to 50 recovers epoch 14 exactly on this run, and so does a naive
  stop-at-first-increase, the lab records in a comment that this is luck on this run rather
  than a property of the rule. Harness check 8 re-runs the capacity and decomposition
  shapes at three data seeds the lesson never quotes, so no directional claim rests on the
  quoted seed.

# Requirements

`requirements.txt` pins the three packages this lab imports directly, at
the exact versions the captured output in `expected-output/` was produced
with:

```
numpy==2.5.2
scikit-learn==1.9.0
pytest==9.1.1
```

Installing scikit-learn also pulls in scipy, joblib and threadpoolctl as
its own dependencies. This lab imports none of them directly and does not
pin them; the versions present during capture are recorded in
`../expected-output/FIELDS.md`.

## Why the versions are pinned exactly

Every dataset here comes from a seeded `numpy.random.default_rng`, and
NumPy's documentation is explicit that `Generator` makes no promise of
stream compatibility between versions. A different NumPy can legitimately
produce a different stream from the same seed, and every measured figure
would move.

What does not depend on the pins: the irreducible floor of 4.0000, which
is arithmetic; the count of 25 features at degree 24; the monotonicity of
training error in capacity and in the penalty; the fact that the three
parts of the decomposition sum to the total; and the *shape* of every
result. Harness check 8 re-runs the shape claims at three data seeds the
lesson never quotes, precisely so the distinction is enforced rather than
asserted.

`../expected-output/FIELDS.md` separates the two categories in full, and
it is worth reading before you conclude that a mismatch is a bug.

## Installing

From the lab directory:

```bash
python3 -m venv .venv
.venv/bin/pip install -r requirements/requirements.txt
```

The install step needs the network. Everything after it is offline: every
dataset in this lab is generated on the spot from a seeded generator, and
nothing is downloaded.

## Free and open-source status

All three packages are free and open source — NumPy and scikit-learn under
the BSD 3-Clause licence, pytest under the MIT licence. There is no paid
tier, no account and no API key anywhere in this lab.

numpy==2.5.2
scikit-learn==1.9.0
pytest==9.1.1

# Day 145 lab brief — Two Ways to Be Wrong

There are exactly two ways a model can be wrong, and they are not two ends
of one dial. They are two different quantities, they respond to completely
different interventions, and this lab measures both directly.

## The claim you are here to measure

> Underfitting is bias. Overfitting is variance. They add up, with
> irreducible noise, to the error you actually observe.

That is usually presented as a picture and left alone, because measuring
its terms needs something you never have: many independent training sets
and knowledge of the true function. In a lab you have both.

Exercise 4 fits **200 models to 200 independent training sets**, predicts
the same fixed grid with all of them, and asks two questions. How far is
the *average* prediction from the truth? That is bias. How far do the
individual predictions scatter from their own average? That is variance.

| degree | bias² | variance | noise | predicted | observed |
| --- | --- | --- | --- | --- | --- |
| 1 | 4.2985 | 0.7112 | 4.0000 | 9.0097 | 9.0295 |
| 3 | 0.0033 | 0.8399 | 4.0000 | 4.8432 | 4.8431 |
| 12 | 2803.5354 | 452183.1336 | 4.0000 | 454990.6691 | 455027.8625 |

The last two columns are the point. **The three parts add up to the error
that was actually observed**, at every capacity, to within one percent.
This is an identity, not an analogy, and you are going to check it.

## The measurement that changes what people do

```text
      n    degree 1     degree 4        degree 24
     15      8.5023      4.9218      215413.2388
   2000      8.2393      3.9880           4.0055
```

A hundred and thirty times more data took the overfit model from 215,413
to 4.0055 — which is exactly the irreducible floor — and the underfit
model from 8.5023 to 8.2393.

**More data cures one failure completely and the other not at all.** If
you are about to spend three months labelling, this is the number that
decides whether it is worth it.

## The finding that is easiest to miss

```text
   deg    train MSE        test MSE             gap
     1      11.3217          9.8274         -1.4942
```

The gap is **negative**. The degree-1 model scores better on data it has
never seen than on the data it was fitted to.

That is not a broken split. A model too rigid to chase the noise in its
training set has no noise-chasing to be flattered by, so its training
score carries none of the usual optimism. It is the signature of
underfitting — and an engineer who sees it usually assumes a leak and
applies more regularisation, which is precisely the wrong direction.

## Two things this lab found by accident, and kept

**Degree 2 is worse than degree 1 on both terms.** It contains every
degree-1 model as a special case and still measures more bias (4.3342
against 4.2985) *and* double the variance. The true function is odd, so a
quadratic term buys nothing and still has to be estimated. Capacity is not
a single dial running from worse to better.

**The degree-24 model is worse at 25 training rows than at 15.** Test
error rises by a factor of three hundred, to sixty-four million, before
falling to the floor. Degree 24 supplies exactly **25 features**, so at 25
rows the system is square: one exact interpolating solution, under no
constraint at all about what happens between the points. At 15 rows there
are more features than rows and least squares returns the minimum-norm
solution, which is quietly a form of regularisation. That peak is the
interpolation threshold, and exercise 3c makes you prove it with a feature
count.

## How to work

1. Build the environment (see the lab `README.md`).
2. Run `.venv/bin/pytest starter -q`. You will see four passes (the
   machinery checks in `test_fitting_lib.py`) and fourteen skips.
3. Replace one `pytest.skip(...)` at a time with real code. The skip text
   names the exact helper and the exact value to assert.
4. Print the measured pair in every exercise. A number you did not print
   is a number you did not look at.
5. When you want the whole measured table at once, run
   `.venv/bin/python3 examples/report_measurements.py`.

Do not run `pytest starter examples` in one invocation. Both directories
define `fitting_lib.py`, `test_fitting_lib.py` and
`test_fitting_claims.py`; pytest aborts on the module-name collision. Run
them separately, always.

## The exercises

| # | What it establishes |
| --- | --- |
| 1 | Training error falls with capacity; test error is U-shaped and explodes |
| 1b | The training curve stops being monotone where the numerics give out |
| 1c | The SIGN of the gap is the diagnostic — negative means underfitting |
| 2 | A penalty rescues the same model class by a factor of 39,588 |
| 2b | The penalty trades training error for test error, and has its own U-curve |
| 3 | More data cures overfitting and does nothing for underfitting |
| 3b | Two good models converge on the irreducible floor from opposite sides |
| 3c | The overfit column peaks exactly where features equal rows |
| 4 | Underfitting is bias; overfitting is variance |
| 4b | The decomposition predicts the error that was observed |
| 4c | A strictly larger model class can be worse on both terms |
| 5 | Training longer improves training error and worsens the model |
| 5b | The generalisation gap grows fivefold while training error falls |
| 5c | The test curve is not a clean U, which is why patience exists |

"""Overfitting and underfitting, decomposed and measured.

The two failures are usually described as a picture: a wiggly line through
every point, a straight line through none. That picture is true and it does
not tell you which one you have, because on your own data you cannot see
the true function.

What you can see is a decomposition. Every squared error splits into three
parts -- bias, variance and irreducible noise -- and the two failures are
one part each:

* **underfitting is bias**: the model class cannot represent the truth, so
  it is wrong in the same direction no matter what data you give it;
* **overfitting is variance**: the model class can represent the truth and
  much else besides, so it chases the noise and lands somewhere different
  for every training set;
* **noise is neither**, and it is the floor nothing can go below.

This module measures all three directly by fitting many models to many
independent training sets, and checks the decomposition against the error
that was actually observed. Everything is deterministic given a seed, and
everything runs on the CPU in seconds.
"""

from __future__ import annotations

import numpy as np

from sklearn.linear_model import LinearRegression, Ridge, SGDRegressor
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler

#: The function the data is actually generated from. A cubic, so a
#: degree-3 model can represent it exactly and a straight line cannot.
NOISE_SD = 2.0


def true_function(x):
    """The relationship the data really has: a cubic, and nothing else."""
    return 0.5 * x**3 - 2.0 * x + 1.0


def irreducible_variance() -> float:
    """The floor. No model of any capacity can score below this."""
    return NOISE_SD**2


def make_data(n: int, seed: int, noise_sd: float = NOISE_SD):
    """n points from the true cubic, plus independent Gaussian noise."""
    rng = np.random.default_rng(seed)
    x = np.sort(rng.uniform(-3.0, 3.0, n))
    y = true_function(x) + rng.normal(0.0, noise_sd, n)
    return x.reshape(-1, 1), y


def polynomial_model(degree: int, alpha: float = 0.0):
    """A polynomial fit of the given degree, optionally ridge-regularised.

    The `StandardScaler` between the features and the estimator is not
    decoration. Raw polynomial features up to degree 24 span many orders
    of magnitude and the normal equations become numerically hopeless;
    without scaling the *training* error starts rising with degree, which
    is a conditioning artefact masquerading as a result.
    """
    estimator = LinearRegression() if alpha == 0.0 else Ridge(alpha=alpha)
    return make_pipeline(PolynomialFeatures(degree), StandardScaler(), estimator)


def mse(model, X, y) -> float:
    """Mean squared error of a fitted model on the given rows."""
    return float(np.mean((model.predict(X) - y) ** 2))


# --------------------------------------------------------------------------
# 1. The capacity sweep
# --------------------------------------------------------------------------


def capacity_sweep(degrees, n_train: int = 25, train_seed: int = 145, test_seed: int = 246):
    """Train and test error at each polynomial degree.

    Rows are ``(degree, train_mse, test_mse, gap)``. Training error falls
    with capacity; test error does not, and where it turns is the whole
    subject.
    """
    X_train, y_train = make_data(n_train, train_seed)
    X_test, y_test = make_data(2000, test_seed)
    rows = []
    for degree in degrees:
        model = polynomial_model(degree).fit(X_train, y_train)
        train = mse(model, X_train, y_train)
        test = mse(model, X_test, y_test)
        rows.append((degree, round(train, 4), round(test, 4), round(test - train, 4)))
    return rows


def best_degree(rows) -> int:
    """The degree with the lowest test error in a capacity sweep."""
    return min(rows, key=lambda row: row[2])[0]


# --------------------------------------------------------------------------
# 2. Regularisation
# --------------------------------------------------------------------------


def regularisation_sweep(alphas, degree: int = 24, n_train: int = 25):
    """Train and test error at each ridge penalty, holding capacity fixed.

    The model class does not change here. Only how much the fit is
    discouraged from using it does -- which is the point: overfitting is
    not a property of the model class alone.
    """
    X_train, y_train = make_data(n_train, 145)
    X_test, y_test = make_data(2000, 246)
    rows = []
    for alpha in alphas:
        model = polynomial_model(degree, alpha=alpha).fit(X_train, y_train)
        rows.append((alpha, round(mse(model, X_train, y_train), 4), round(mse(model, X_test, y_test), 4)))
    return rows


def best_alpha(rows) -> float:
    """The penalty with the lowest test error."""
    return min(rows, key=lambda row: row[2])[0]


# --------------------------------------------------------------------------
# 3. What more data fixes, and what it does not
# --------------------------------------------------------------------------


def data_sweep(sizes, degrees=(1, 4, 24), test_seed: int = 246):
    """Test error against training-set size, at three capacities.

    Rows are ``(n, {degree: test_mse})``. The three columns behave
    completely differently, and the difference is the practical content of
    the whole lesson.
    """
    X_test, y_test = make_data(2000, test_seed)
    rows = []
    for n in sizes:
        X_train, y_train = make_data(n, 900 + n)
        scores = {}
        for degree in degrees:
            model = polynomial_model(degree).fit(X_train, y_train)
            scores[degree] = round(mse(model, X_test, y_test), 4)
        rows.append((n, scores))
    return rows


# --------------------------------------------------------------------------
# 4. The decomposition itself
# --------------------------------------------------------------------------


def bias_variance(degree: int, n_train: int = 25, datasets: int = 200, grid: int = 200):
    """Measure bias squared, variance and noise for one model class.

    Fit ``datasets`` models, each to its own independent training set, and
    predict the same fixed grid of query points with all of them. Then:

    * **bias squared** is how far the *average* prediction sits from the
      truth -- an error the model class makes every time;
    * **variance** is how much the predictions scatter around their own
      average -- an error that changes with the training set;
    * **noise** is the irreducible term, known here because we generated it.

    Returns a dict with all three, their sum, and the squared error
    actually observed against freshly noisy targets. Those last two should
    agree, and checking that they do is what makes this a measurement
    rather than a recitation.
    """
    query = np.linspace(-3.0, 3.0, grid).reshape(-1, 1)
    truth = true_function(query.ravel())

    predictions = np.empty((datasets, grid))
    for i in range(datasets):
        X_train, y_train = make_data(n_train, 10_000 + i)
        model = polynomial_model(degree).fit(X_train, y_train)
        predictions[i] = model.predict(query)

    mean_prediction = predictions.mean(axis=0)
    bias_squared = float(np.mean((mean_prediction - truth) ** 2))
    variance = float(np.mean(predictions.var(axis=0)))
    noise = irreducible_variance()

    observed_rng = np.random.default_rng(7)
    noisy_targets = truth + observed_rng.normal(0.0, NOISE_SD, predictions.shape)
    observed = float(np.mean((predictions - noisy_targets) ** 2))

    return {
        "bias_squared": round(bias_squared, 4),
        "variance": round(variance, 4),
        "noise": round(noise, 4),
        "predicted_total": round(bias_squared + variance + noise, 4),
        "observed": round(observed, 4),
    }


def decomposition_table(degrees, n_train: int = 25, datasets: int = 200):
    """The decomposition at several capacities, as one table."""
    return [(degree, bias_variance(degree, n_train, datasets)) for degree in degrees]


# --------------------------------------------------------------------------
# 5. Early stopping
# --------------------------------------------------------------------------


def training_history(epochs: int = 600, degree: int = 14, n_train: int = 25, seed: int = 0):
    """Train and test error after every epoch of gradient descent.

    Capacity is fixed and the data is fixed. The only thing that changes is
    how long the fit has been allowed to run -- which turns out to be a
    capacity knob of its own.
    """
    X_train, y_train = make_data(n_train, 145)
    X_test, y_test = make_data(2000, 246)

    features = PolynomialFeatures(degree)
    scaler = StandardScaler()
    P_train = scaler.fit_transform(features.fit_transform(X_train))
    P_test = scaler.transform(features.transform(X_test))

    model = SGDRegressor(
        learning_rate="constant", eta0=0.003, penalty=None, random_state=seed
    )
    train_history, test_history = [], []
    for _epoch in range(epochs):
        model.partial_fit(P_train, y_train)
        train_history.append(float(np.mean((model.predict(P_train) - y_train) ** 2)))
        test_history.append(float(np.mean((model.predict(P_test) - y_test) ** 2)))
    return train_history, test_history


def is_monotonically_decreasing(values, tolerance: float = 1e-9) -> bool:
    """Whether a sequence never increases, within a numerical tolerance."""
    return all(a >= b - tolerance for a, b in zip(values, values[1:]))


def first_increase(values) -> int:
    """The index of the first value larger than the one before it."""
    for i, (a, b) in enumerate(zip(values, values[1:]), start=1):
        if b > a:
            return i
    return len(values)


def stop_with_patience(values, patience: int) -> int:
    """Index of the epoch a patience-based early stop would have chosen.

    Stops once ``patience`` consecutive epochs have failed to improve on
    the best seen so far, and returns the best epoch rather than the epoch
    it stopped at -- which is what a real implementation restores.
    """
    best_index, best_value, waited = 0, values[0], 0
    for i, value in enumerate(values[1:], start=1):
        if value < best_value:
            best_index, best_value, waited = i, value, 0
        else:
            waited += 1
            if waited >= patience:
                break
    return best_index

"""Fourteen exercises in the two ways a model can be wrong.

Read `00_brief.md` first. Each function below is a `pytest.skip` naming
exactly what to build and what to assert; replace the skip with real code.
`fitting_lib.py` is complete -- it is the machinery, not the exercise.

Run this suite on its own:

    .venv/bin/pytest starter -q

Never run `pytest starter examples` in one invocation: both directories
define modules with the same names and pytest aborts on the collision.
"""

import numpy as np  # noqa: F401  (you will need it)
import pytest

import fitting_lib as f  # noqa: F401  (you will need it)


def test_01_training_error_falls_with_capacity_and_test_error_does_not():
    pytest.skip(
        "Assert f.capacity_sweep([1, 2, 3, 4, 6, 8, 10, 14, 18, 24]) equals "
        "the ten rows in expected-output/measured-values.txt, from (1, "
        "11.3217, 9.8274, -1.4942) to (24, 1.0321, 226667.4689, "
        "226666.4368). Assert f.best_degree is 4 and that the test error at "
        "degree 24 is more than 40000 times the test error at degree 4."
    )


def test_01b_training_error_falls_until_the_numerics_give_out():
    pytest.skip(
        "Take the training column of that sweep. Assert it is monotonically "
        "decreasing through degree 14 but NOT over the whole range, and that "
        "the wobble after degree 14 is 0.0636. That wobble is numerical, not "
        "statistical: 25 training rows, and degree 24 supplies exactly 25 "
        "features. Say so in a comment."
    )


def test_01c_the_generalisation_gap_is_the_diagnostic():
    pytest.skip(
        "Take the gap column. Assert the gap is NEGATIVE at degrees 1 and 2 "
        "-- test error below training error, because a model too rigid to "
        "chase noise has none to be flattered by. Assert it is between 0 and "
        "3 at the best degree, and above 200000 at degree 24. The gap tells "
        "you which failure you have without ever seeing the true function."
    )


def test_02_a_penalty_rescues_the_same_model_class():
    pytest.skip(
        "Assert f.regularisation_sweep([0.0, 1e-6, 1e-3, 0.1, 1.0, 10.0, "
        "100.0]) matches the captured rows, that f.best_alpha is 1.0, and "
        "that the ratio of the unpenalised test error to the best one is "
        "39588. The degree never changed. Only how strongly the fit was "
        "discouraged from using it."
    )


def test_02b_the_penalty_trades_training_error_for_test_error():
    pytest.skip(
        "Assert the training column rises monotonically with the penalty -- "
        "always true, since the penalty can only make the training fit "
        "worse. Assert the test column is U-shaped too, with alpha 1.0 "
        "beating 10.0 beating 100.0. Then assert the best-regularised "
        "degree-24 model lands within 0.3 of the best unregularised "
        "degree-4 model from f.capacity_sweep([4])."
    )


def test_03_more_data_cures_overfitting_and_does_nothing_for_underfitting():
    pytest.skip(
        "Assert f.data_sweep([15, 25, 50, 100, 400, 2000]) matches the "
        "captured rows. Then assert the two things that matter: the "
        "degree-1 column spans only 0.6227 across a 133-fold increase in "
        "data and starts and ends within 0.3 of itself, while the degree-24 "
        "column falls by more than seven orders of magnitude."
    )


def test_03b_both_good_models_converge_to_the_irreducible_floor():
    pytest.skip(
        "Assert f.irreducible_variance() is 4.0. At n=2000, assert degree 4 "
        "is within 0.02 of that floor and degree 24 within 0.01 -- they "
        "arrive from opposite sides, one from below the floor's own noise "
        "and one from far above. Assert the degree-1 model is still more "
        "than twice the floor at n=2000, and will be at any n."
    )


def test_03c_the_overfit_column_peaks_at_the_interpolation_threshold():
    pytest.skip(
        "The degree-24 column is NOT monotone in n: it is worse at n=25 "
        "than at n=15. Assert that, and assert n=25 is its maximum. Then "
        "explain it: use sklearn.preprocessing.PolynomialFeatures(24) on a "
        "dummy array to assert it produces exactly 25 features, and assert "
        "the peak row's n equals that count. Features equal to rows is the "
        "worst-conditioned case there is."
    )


def test_04_underfitting_is_bias_and_overfitting_is_variance():
    pytest.skip(
        "Call f.bias_variance at degrees 1, 3 and 12. Assert bias squared "
        "is 4.2985, 0.0033 and 2803.5354 and variance is 0.7112, 0.8399 and "
        "452183.1336. Then assert the shape of the story: at degree 1 bias "
        "exceeds variance sixfold, at degree 3 the bias is below 0.01 "
        "because the true function IS a cubic, and at degree 12 variance "
        "exceeds bias a hundredfold."
    )


def test_04b_the_decomposition_predicts_the_error_that_was_observed():
    pytest.skip(
        "For each degree in (1, 2, 3, 4, 6, 8, 12), assert the predicted "
        "total agrees with the observed squared error to within 1.1 percent, "
        "and that the worst disagreement across all seven is exactly "
        "0.01003. Also assert the three parts sum to the total to within "
        "0.0002 -- not exactly, because each part is stored already rounded "
        "to four places."
    )


def test_04c_a_bigger_model_class_is_not_always_less_biased():
    pytest.skip(
        "Assert degree 2 has bias squared 4.3342 against degree 1's 4.2985 "
        "-- MORE bias from a strictly larger model class. Assert it also "
        "carries roughly double the variance, and that its predicted total "
        "is therefore worse on both counts. The true function is odd, so a "
        "quadratic term can buy nothing and still costs. Capacity is not a "
        "single dial."
    )


def test_05_training_longer_makes_training_error_better_and_the_model_worse():
    pytest.skip(
        "Call f.training_history(). Assert both histories have 600 entries, "
        "that the training history is monotonically decreasing over all 600 "
        "epochs, and that it runs 7.3906 -> 2.4744. Then assert the test "
        "error bottoms at index 13 (epoch 14) with 5.4555 and is 5.8978 at "
        "epoch 600. Nothing about the model or the data changed; only how "
        "long the fit ran."
    )


def test_05b_the_generalisation_gap_grows_while_training_error_falls():
    pytest.skip(
        "Assert the gap between test and training error is 0.6771 at epoch "
        "1 and 3.4234 at epoch 600, a factor of more than five. Training "
        "time is a capacity knob, and the gap is what it is spending."
    )


def test_05c_the_test_curve_is_not_a_clean_u_which_is_why_patience_exists():
    pytest.skip(
        "After the best epoch, assert the test error rises to 7.1435 and "
        "then partly recovers to 5.8978 -- without ever again beating the "
        "5.4555 it reached at epoch 14, so assert min(after) is still above "
        "it and that 599 of the 600 epochs are worse than the best. Then "
        "assert every patience in (5, 10, 20, 50) recovers epoch 14 exactly, "
        "and that f.first_increase does too. Note in a comment that this is "
        "luck on this run rather than a property of the rule -- a curve that "
        "wanders like this one can defeat a naive stop-on-first-rise."
    )

"""Machinery checks: the helpers behave, before any claim is made.

These four tests are solved in both `starter/` and `examples/`. They exist
so that a broken helper reports itself as a broken helper rather than as a
surprising scientific result.
"""

import numpy as np
import pytest

import fitting_lib as f


def test_the_data_generator_is_the_true_function_plus_noise():
    X, y = f.make_data(2000, seed=1)
    residuals = y - f.true_function(X.ravel())
    # The noise is centred and has the standard deviation it claims.
    assert abs(float(residuals.mean())) < 0.15
    assert abs(float(residuals.std()) - f.NOISE_SD) < 0.1
    assert f.irreducible_variance() == f.NOISE_SD**2 == 4.0
    # Two calls at one seed agree; two seeds do not.
    assert np.array_equal(f.make_data(50, 3)[1], f.make_data(50, 3)[1])
    assert not np.array_equal(f.make_data(50, 3)[1], f.make_data(50, 4)[1])


def test_a_degree_three_model_can_represent_the_truth_exactly():
    # With no noise at all, a cubic fit should be essentially perfect and
    # a straight line should not.
    X, y = f.make_data(200, seed=2, noise_sd=0.0)
    cubic = f.polynomial_model(3).fit(X, y)
    line = f.polynomial_model(1).fit(X, y)
    assert f.mse(cubic, X, y) < 1e-12
    # A straight line cannot, and is left with 4.174 of pure bias.
    assert round(f.mse(line, X, y), 3) == 4.174


def test_scaling_is_what_keeps_the_high_degree_fit_numerically_sane():
    from sklearn.linear_model import LinearRegression
    from sklearn.pipeline import make_pipeline
    from sklearn.preprocessing import PolynomialFeatures

    X, y = f.make_data(25, seed=145)
    scaled = f.polynomial_model(24).fit(X, y)
    unscaled = make_pipeline(PolynomialFeatures(24), LinearRegression()).fit(X, y)
    # Both fit the training data; the scaled pipeline fits it better,
    # because the unscaled normal equations are badly conditioned.
    assert f.mse(scaled, X, y) < f.mse(unscaled, X, y)


def test_the_stopping_helpers_agree_on_a_hand_checkable_sequence():
    values = [10.0, 8.0, 5.0, 6.0, 7.0, 9.0]
    assert f.first_increase(values) == 3
    assert not f.is_monotonically_decreasing(values)
    assert f.is_monotonically_decreasing([5.0, 4.0, 4.0, 3.0])
    # Patience 2 waits two non-improving epochs, then returns the best.
    assert f.stop_with_patience(values, 2) == 2
    # A sequence that keeps improving is never stopped early.
    assert f.stop_with_patience([5.0, 4.0, 3.0, 2.0], 2) == 3

#!/usr/bin/env bash
# Day 145 lab harness: "Two Ways to Be Wrong"
#
# Prints "N checks, M failure(s)" and exits 0 only when M is zero.
set -u

LAB_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)"
cd "$LAB_DIR"

PYTHON="${PYTHON:-.venv/bin/python3}"
PYTEST="${PYTEST:-.venv/bin/pytest}"

find . -path ./.venv -prune -o -type d -name '__pycache__' -exec rm -rf -- {} + 2>/dev/null
rm -rf .pytest_cache

CHECKS=0
FAILURES=0
ok() { CHECKS=$((CHECKS + 1)); echo "  ok: $1"; }
fail() { CHECKS=$((CHECKS + 1)); FAILURES=$((FAILURES + 1)); echo "  FAIL: $1"; }

if [ ! -x "$PYTHON" ]; then
  echo "No lab .venv found at $PYTHON."
  echo "Run: python3 -m venv .venv && .venv/bin/pip install -r requirements/requirements.txt"
  exit 2
fi

echo "1. Installed versions match requirements/requirements.txt"
VERSION_CHECK=$("$PYTHON" - <<'PYEOF'
import numpy, sklearn, pytest
print("numpy", numpy.__version__)
print("scikit-learn", sklearn.__version__)
print("pytest", pytest.__version__)
PYEOF
)
echo "$VERSION_CHECK" | sed 's/^/    /'
while read -r pkg pin; do
  pin_version="${pin#*==}"
  installed=$(echo "$VERSION_CHECK" | awk -v p="$pkg" '$1==p {print $2}')
  if [ "$installed" = "$pin_version" ]; then ok "$pkg $installed matches the pin"
  else fail "$pkg installed=$installed pinned=$pin_version"; fi
done < <(sed 's/==/ ==/' requirements/requirements.txt)

echo ""
echo "2. Every published claim, reproduced directly (no pytest involved)"
DIRECT_CHECK=$("$PYTHON" - <<'PYEOF'
import sys
sys.path.insert(0, "examples")
import numpy as np
import fitting_lib as f

errors = []
def expect(label, got, want):
    if got != want:
        errors.append(f"{label}: expected {want}, got {got}")

# 1. Capacity
rows = f.capacity_sweep([1, 2, 3, 4, 6, 8, 10, 14, 18, 24])
expect("capacity sweep", rows, [
    (1, 11.3217, 9.8274, -1.4942), (2, 11.3173, 9.8801, -1.4372),
    (3, 2.7076, 6.123, 3.4154), (4, 2.4964, 5.4911, 2.9948),
    (6, 1.9569, 15.8217, 13.8648), (8, 1.701, 26.1708, 24.4697),
    (10, 1.357, 528.4798, 527.1227), (14, 0.9685, 31307.2782, 31306.3097),
    (18, 1.0037, 75539.3618, 75538.3581), (24, 1.0321, 226667.4689, 226666.4368)])
expect("best degree", f.best_degree(rows), 4)
train = [r[1] for r in rows]
expect("train monotone through degree 14", f.is_monotonically_decreasing(train[:8]), True)
expect("train monotone overall", f.is_monotonically_decreasing(train), False)
expect("post-14 wobble", round(max(train[7:]) - min(train[7:]), 4), 0.0636)
gaps = {d: g for d, _t, _e, g in rows}
if not (gaps[1] < 0 and gaps[2] < 0):
    errors.append("the underfitting gap was not negative")
if not (0 < gaps[4] < 3.0):
    errors.append("the best model's gap was not small and positive")

# 2. Regularisation
reg = f.regularisation_sweep([0.0, 1e-6, 1e-3, 0.1, 1.0, 10.0, 100.0])
expect("regularisation sweep", reg, [
    (0.0, 1.0321, 226667.4689), (1e-06, 1.2031, 130776.6548),
    (0.001, 1.5741, 128.3127), (0.1, 2.2259, 15.8339),
    (1.0, 2.7461, 5.7257), (10.0, 3.9689, 6.1559), (100.0, 6.28, 6.784)])
expect("best alpha", f.best_alpha(reg), 1.0)
expect("rescue factor", round(reg[0][2] / min(r[2] for r in reg), 0), 39588.0)
reg_train = [r[1] for r in reg]
if not all(a < b for a, b in zip(reg_train, reg_train[1:])):
    errors.append("training error did not rise monotonically with the penalty")

# 3. Data
data = f.data_sweep([15, 25, 50, 100, 400, 2000])
expect("data sweep", data, [
    (15, {1: 8.5023, 4: 4.9218, 24: 215413.2388}),
    (25, {1: 8.862, 4: 6.1904, 24: 64631547.2994}),
    (50, {1: 8.2457, 4: 4.2661, 24: 6070.3302}),
    (100, {1: 8.3583, 4: 4.2934, 24: 5.3571}),
    (400, {1: 8.3007, 4: 3.9958, 24: 4.3139}),
    (2000, {1: 8.2393, 4: 3.988, 24: 4.0055})])
expect("irreducible floor", f.irreducible_variance(), 4.0)
underfit = [s[1] for _n, s in data]
overfit = [s[24] for _n, s in data]
expect("underfit range", round(max(underfit) - min(underfit), 4), 0.6227)
if overfit[1] / overfit[-1] < 1e7:
    errors.append("more data did not cure the overfit model by seven orders of magnitude")
if abs(data[-1][1][24] - 4.0) > 0.01 or abs(data[-1][1][4] - 4.0) > 0.02:
    errors.append("the good models did not converge to the irreducible floor")
if overfit[1] != max(overfit):
    errors.append("the overfit column did not peak at the interpolation threshold")

# 4. Decomposition
worst = 0.0
for degree in (1, 2, 3, 4, 6, 8, 12):
    r = f.bias_variance(degree)
    rel = abs(r["predicted_total"] - r["observed"]) / r["observed"]
    worst = max(worst, rel)
    if rel >= 0.011:
        errors.append(f"degree {degree}: decomposition off by {rel:.4f}")
    if abs(r["bias_squared"] + r["variance"] + r["noise"] - r["predicted_total"]) > 0.0002:
        errors.append(f"degree {degree}: the parts did not sum to the total")
expect("worst decomposition disagreement", round(worst, 5), 0.01003)
one, three, twelve = f.bias_variance(1), f.bias_variance(3), f.bias_variance(12)
expect("degree 1 bias squared", one["bias_squared"], 4.2985)
expect("degree 1 variance", one["variance"], 0.7112)
expect("degree 3 bias squared", three["bias_squared"], 0.0033)
expect("degree 12 variance", twelve["variance"], 452183.1336)
if one["bias_squared"] / one["variance"] <= 6:
    errors.append("bias did not dominate at degree 1")
if twelve["variance"] / twelve["bias_squared"] <= 100:
    errors.append("variance did not dominate at degree 12")
two = f.bias_variance(2)
expect("degree 2 bias squared", two["bias_squared"], 4.3342)
if two["bias_squared"] <= one["bias_squared"]:
    errors.append("degree 2 did not have more bias than degree 1")

# 5. Early stopping
train_hist, test_hist = f.training_history()
expect("epochs", len(train_hist), 600)
expect("training monotone", f.is_monotonically_decreasing(train_hist), True)
expect("first training error", round(train_hist[0], 4), 7.3906)
expect("last training error", round(train_hist[-1], 4), 2.4744)
best = int(np.argmin(test_hist))
expect("best epoch index", best, 13)
expect("best test error", round(test_hist[best], 4), 5.4555)
expect("test error at 600", round(test_hist[-1], 4), 5.8978)
expect("worst after best", round(max(test_hist[best + 1:]), 4), 7.1435)
expect("epochs worse than best", sum(1 for v in test_hist if v > min(test_hist)), 599)
expect("first gap", round(test_hist[0] - train_hist[0], 4), 0.6771)
expect("last gap", round(test_hist[-1] - train_hist[-1], 4), 3.4234)
for patience in (5, 10, 20, 50):
    expect(f"patience {patience}", f.stop_with_patience(test_hist, patience), best)

if errors:
    for e in errors:
        print("ERROR:", e)
    sys.exit(1)
print("all direct checks passed")
PYEOF
)
if echo "$DIRECT_CHECK" | grep -q "all direct checks passed"; then
  ok "exercises 1-5 reproduced directly against fitting_lib, no pytest involved"
else
  fail "direct library checks failed"
  echo "$DIRECT_CHECK" | sed 's/^/    /'
fi

echo ""
echo "3. examples/ passes in full"
EXAMPLES_OUT=$("$PYTEST" examples -q 2>&1)
if echo "$EXAMPLES_OUT" | tail -1 | grep -qE "^18 passed"; then
  ok "pytest examples -q -> 18 passed"
else
  fail "pytest examples -q did not report 18 passed"
  echo "$EXAMPLES_OUT" | tail -20 | sed 's/^/    /'
fi

echo ""
echo "4. starter/ is an untouched skeleton"
STARTER_OUT=$("$PYTEST" starter -q 2>&1)
if echo "$STARTER_OUT" | tail -1 | grep -qE "4 passed, 14 skipped"; then
  ok "pytest starter -q -> 4 passed, 14 skipped (the machinery checks pass; the fourteen exercises are stubs)"
else
  fail "pytest starter -q did not report 4 passed, 14 skipped"
  echo "$STARTER_OUT" | tail -20 | sed 's/^/    /'
fi

echo ""
echo "5. pytest examples starter (one invocation) aborts on the module-name collision"
COMBINED_OUT=$("$PYTEST" examples starter 2>&1)
if echo "$COMBINED_OUT" | grep -q "import file mismatch"; then
  ok "combined invocation reports import file mismatch, as documented -- never run starter and examples together"
else
  fail "combined invocation did not fail with import file mismatch as expected"
fi

echo ""
echo "6. The report reproduces the captured table exactly"
REPORT_OUT=$("$PYTHON" examples/report_measurements.py 2>&1)
if [ "$REPORT_OUT" = "$(cat expected-output/measured-values.txt)" ]; then
  ok "report_measurements.py output is byte-identical to expected-output/measured-values.txt"
else
  fail "report_measurements.py drifted from expected-output/measured-values.txt"
  echo "$REPORT_OUT" | diff - expected-output/measured-values.txt | head -20 | sed 's/^/    /'
fi

echo ""
echo "7. Proof the harness can fail"
SCRATCH=$(mktemp -d "${TMPDIR:-/tmp}/d145-scratch.XXXXXX")
cp examples/*.py "$SCRATCH"/
SCRATCH_OUT=$("$PYTEST" "$SCRATCH" -q 2>&1)
if echo "$SCRATCH_OUT" | tail -1 | grep -qE "^18 passed"; then
  ok "scratch copy of examples/ passes before it is broken"
else
  fail "scratch copy did not pass before being broken: $(echo "$SCRATCH_OUT" | tail -3)"
fi
"$PYTHON" - "$SCRATCH/test_fitting_claims.py" <<'PYEOF'
import sys
path = sys.argv[1]
text = open(path).read()
needle = "assert f.best_degree(rows) == 4"
replacement = "assert f.best_degree(rows) == 24"
assert needle in text, "could not find the assertion to break"
open(path, "w").write(text.replace(needle, replacement, 1))
PYEOF
BROKEN_OUT=$("$PYTEST" "$SCRATCH" -q 2>&1)
BROKEN_STATUS=$?
if [ "$BROKEN_STATUS" -ne 0 ] && echo "$BROKEN_OUT" | grep -q "test_01_training_error_falls_with_capacity_and_test_error_does_not"; then
  ok "breaking exercise 1's assertion produces a non-zero exit and names the failing test"
else
  fail "broken copy did not fail as expected (exit=$BROKEN_STATUS)"
fi
rm -rf "$SCRATCH"

echo ""
echo "8. The shape of every result survives a different data seed"
SHAPE=$("$PYTHON" - <<'PYEOF'
import sys
sys.path.insert(0, "examples")
import numpy as np
import fitting_lib as f

problems = []
for seed in (301, 302, 303):
    rows = f.capacity_sweep([1, 3, 4, 8, 24], train_seed=seed, test_seed=seed + 50)
    test = {d: t for d, _tr, t, _g in rows}
    train = {d: tr for d, tr, _t, _g in rows}
    if train[24] >= train[1]:
        problems.append(f"seed {seed}: training error did not fall with capacity")
    if test[24] <= test[4]:
        problems.append(f"seed {seed}: the degree-24 model did not overfit")
    if test[1] <= test[4]:
        problems.append(f"seed {seed}: the degree-1 model did not underfit")

# And the decomposition's shape: bias dominant when rigid, variance when not.
low = f.bias_variance(1, datasets=60)
high = f.bias_variance(12, datasets=60)
if low["bias_squared"] <= low["variance"]:
    problems.append("bias did not dominate for the rigid model at a smaller sample")
if high["variance"] <= high["bias_squared"]:
    problems.append("variance did not dominate for the flexible model at a smaller sample")

if problems:
    for p in problems:
        print("ERROR:", p)
else:
    print("every shape held")
PYEOF
)
if [ "$SHAPE" = "every shape held" ]; then
  ok "underfitting, overfitting and the decomposition's shape hold at seeds the lesson does not quote"
else
  fail "a shape failed beyond the quoted seed"
  echo "$SHAPE" | sed 's/^/    /'
fi

echo ""
echo "9. Offline, and nothing left behind"
if ! grep -rInE "https?://" examples/*.py starter/*.py > /dev/null 2>&1; then
  ok "no URLs inside examples/ or starter/ source -- this lab reaches no network"
else
  fail "found a URL inside examples/ or starter/"
fi
if [ -z "$(find . -path ./.venv -prune -o -type d -name '__pycache__' -print 2>/dev/null)" ]; then
  ok "no __pycache__ left behind"
else
  find . -path ./.venv -prune -o -type d -name '__pycache__' -exec rm -rf -- {} + 2>/dev/null
  ok "no __pycache__ left behind (cleaned during this run)"
fi
if [ ! -d .pytest_cache ]; then
  ok "no .pytest_cache left behind"
else
  rm -rf .pytest_cache
  ok "no .pytest_cache left behind (cleaned during this run)"
fi

echo ""
echo "---------------------------------------------------------------"
echo "$CHECKS checks, $FAILURES failure(s)"
if [ "$FAILURES" -ne 0 ]; then exit 1; fi
exit 0