Hands-on lab — Day 144: Train, Validation, and Test Splits

#!/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
from sklearn.linear_model import LogisticRegression  # noqa: E402

import splits_lib as s  # noqa: E402


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


def main() -> None:
    print("Day 144 -- train, validation and test splits, measured")
    print("=" * 53)

    rule("1. Why three sets: selecting on a set is fitting to it")
    standard_error = s.proportion_standard_error(0.5, 500)
    print("  K candidates, each a coin flip, 500-row validation and 500-row test")
    print(f"  averaged over 400 replications; standard error of an accuracy is {standard_error:.4f}")
    print("     K   best-val   its-test   optimism   in SEs")
    rows = s.selection_bias_curve([1, 2, 5, 10, 25, 50, 100, 500, 1000])
    for k, validation, test, optimism in rows:
        print(f"  {k:5d}    {validation:.4f}     {test:.4f}    {optimism:+.4f}    {optimism / standard_error:5.2f}")
    print("  the test column never moves: it was never selected on")

    rule("1b. The optimism is the expected maximum of K noise draws")
    print("     K   measured   E max of K normals   sqrt(2 ln K)")
    for k, _v, _t, optimism in rows[1:]:
        print(
            f"  {k:5d}     {optimism / standard_error:5.2f}              "
            f"{s.expected_max_of_normals(k):5.2f}          {s.sqrt_two_log_k(k):5.2f}"
        )
    print("  the closed-form approximation sits above the truth at every K here")

    rule("2. Stratification: a rare class and a small test set")
    X, y = s.rare_class_dataset()
    random_rates, stratified_rates, empty = s.split_positive_rates(X, y)
    print(f"  population positive rate : {float(y.mean()):.4f}")
    print(f"  random split     : {s.spread(random_rates)}")
    print(f"  stratified split : {s.spread(stratified_rates)}")
    print(f"  random splits whose test half held ZERO positives: {empty} of 500")

    rule("3. Groups: when the row is not the unit")
    Xg, yg, groups = s.grouped_dataset()
    rowwise, group_aware = s.rowwise_vs_group_split(Xg, yg, groups)
    print("  50 people, 20 rows each; each person's label is a coin flip")
    print(f"  row-wise random split : {rowwise:.4f}")
    print(f"  group-aware split     : {group_aware:.4f}")
    print(f"  accuracy invented     : {rowwise - group_aware:+.4f}")
    print(f"  people appearing in BOTH halves of a row-wise split: {s.groups_shared_between_halves(groups)} of 50")

    rule("4. Time: when the data has a direction")
    temporal = s.temporal_inflation_over_constructions()
    inflation = [gap for _s, _sh, _ch, _b, gap in temporal]
    print("  20 independently generated series, six regimes each")
    print(f"  shuffled split, mean       : {np.mean([r[1] for r in temporal]):.4f}")
    print(f"  chronological split, mean  : {np.mean([r[2] for r in temporal]):.4f}")
    print(f"  majority baseline, mean    : {np.mean([r[3] for r in temporal]):.4f}")
    print(f"  inflation: mean {np.mean(inflation):.4f}  sd {np.std(inflation):.4f}  "
          f"min {min(inflation):+.4f}  max {max(inflation):+.4f}")
    print(f"  constructions where shuffling won: {sum(1 for g in inflation if g > 0)} of 20")
    print("  the direction is universal; the size varies by a factor of sixteen")

    rule("5. One holdout, or many folds")
    Xw, yw = s.weak_signal_dataset()
    holdout, cross = s.holdout_vs_cross_validation(Xw, yw)
    print("  same data, same model; only which rows landed where changes")
    print(f"  single holdout : {s.spread(holdout)}")
    print(f"  5-fold         : {s.spread(cross)}")
    print(f"  holdout swings {max(holdout) - min(holdout):.4f} across 200 seeds")
    print(f"  cross-validation is {np.std(holdout) / np.std(cross):.4f} times steadier")

    rule("6. How big does the test set need to be?")
    print("       n   theory SE   measured sd   95 percent half-width")
    for n, theory, measured, half in s.test_size_table([50, 100, 200, 500, 1000, 5000]):
        print(f"  {n:6d}      {theory:.4f}        {measured:.4f}                 +/-{half:.4f}")
    print(f"  rows needed for +/-0.02 at an accuracy of 0.85 : {s.rows_needed_for_precision(0.85, 0.02)}")
    print(f"  rows needed for +/-0.01 at an accuracy of 0.85 : {s.rows_needed_for_precision(0.85, 0.01)}")

    rule("7. The rule, made mechanical")
    model = LogisticRegression(max_iter=1000).fit(Xw, yw)
    gate = s.GatedTestSet(Xw, yw)
    print(f"  first evaluation  : {gate.evaluate(model):.4f}")
    try:
        gate.evaluate(model)
        print("  second evaluation : NO ERROR RAISED")
    except s.TestSetTouchedTwice as exc:
        print(f"  second evaluation : {type(exc).__name__}")
        print(f"    {exc}")


if __name__ == "__main__":
    main()

"""Splitting, measured: what each way of cutting a dataset actually buys.

Three sets, not two, and the reason is arithmetic rather than convention.
A validation set you select on is a set you have fitted to, and this
module measures the resulting optimism directly -- it turns out to be
exactly the expected maximum of K noise draws, which is a quantity you can
compute.

The rest measures the four ways a split goes wrong: not stratifying when
the class is rare, splitting rows when the unit is a person, splitting
randomly when the data has a direction in time, and reading a trend off
one holdout when a holdout's own spread is wider than the trend.

Everything here is deterministic given a seed.
"""

from __future__ import annotations

import numpy as np

from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import (
    GroupShuffleSplit,
    StratifiedKFold,
    StratifiedShuffleSplit,
    cross_val_score,
    train_test_split,
)
from sklearn.neighbors import KNeighborsClassifier


def accuracy(y_true, y_pred) -> float:
    return float(np.mean(np.asarray(y_true) == np.asarray(y_pred)))


# --------------------------------------------------------------------------
# 1. Why three sets: selecting on a set is fitting to it
# --------------------------------------------------------------------------


def selection_replicate(k_candidates: int, n: int = 500, seed: int = 0):
    """One replication of "pick the best of K, then check it on a fresh set".

    Each candidate is a fixed prediction vector over 2n rows, split into an
    n-row validation set and an n-row test set. Every candidate has exactly
    zero skill by construction, so any validation score above 0.5 is noise
    and any *selected* validation score above 0.5 is noise you chose.

    Returns ``(best_validation_score, that_candidate_s_test_score)``.
    """
    rng = np.random.default_rng(seed)
    y = rng.integers(0, 2, size=2 * n)
    predictions = rng.integers(0, 2, size=(k_candidates, 2 * n))
    correct = predictions == y
    validation = correct[:, :n].mean(axis=1)
    test = correct[:, n:].mean(axis=1)
    winner = int(np.argmax(validation))
    return float(validation[winner]), float(test[winner])


def selection_bias_curve(k_values, replications: int = 400, n: int = 500):
    """Mean selected-validation score and its test score, for each K.

    Returns rows of ``(k, mean_validation, mean_test, optimism)``. The test
    column is the control: it must stay at chance for every K, because the
    test set was never selected on.
    """
    rows = []
    for k in k_values:
        pairs = [selection_replicate(k, n=n, seed=r) for r in range(replications)]
        validation = float(np.mean([v for v, _t in pairs]))
        test = float(np.mean([t for _v, t in pairs]))
        rows.append((k, round(validation, 4), round(test, 4), round(validation - test, 4)))
    return rows


def proportion_standard_error(p: float, n: int) -> float:
    """The standard error of an accuracy estimated on n rows."""
    return float(np.sqrt(p * (1.0 - p) / n))


def expected_max_of_normals(k: int, draws: int = 20000, seed: int = 7) -> float:
    """E of the maximum of k standard normals, by simulation.

    This is the quantity the selection optimism should equal, once the
    optimism is expressed in standard errors. Simulated rather than
    approximated, because the usual closed form overestimates it -- which
    exercise 1c measures.
    """
    rng = np.random.default_rng(seed)
    return float(np.mean(np.max(rng.standard_normal((draws, k)), axis=1)))


def sqrt_two_log_k(k: int) -> float:
    """The textbook asymptotic for the expected maximum of k normals."""
    return 0.0 if k <= 1 else float(np.sqrt(2.0 * np.log(k)))


# --------------------------------------------------------------------------
# 2. Stratification: a rare class and a small test set
# --------------------------------------------------------------------------


def rare_class_dataset(n: int = 200, rate: float = 0.05, seed: int = 144):
    """A dataset with a genuinely rare positive class."""
    rng = np.random.default_rng(seed)
    y = np.zeros(n, dtype=int)
    y[: int(n * rate)] = 1
    X = rng.normal(size=(n, 3))
    return X, y


def split_positive_rates(X, y, splits: int = 500, test_size: float = 0.25):
    """Positive rate in the test half, over many random and stratified splits.

    Returns ``(random_rates, stratified_rates, random_splits_with_no_positives)``.
    """
    random_rates = []
    stratified_rates = []
    empty = 0
    for seed in range(splits):
        _x_tr, _x_te, _y_tr, y_te = train_test_split(X, y, test_size=test_size, random_state=seed)
        random_rates.append(float(y_te.mean()))
        if y_te.sum() == 0:
            empty += 1
        splitter = StratifiedShuffleSplit(n_splits=1, test_size=test_size, random_state=seed)
        _train, test = next(splitter.split(X, y))
        stratified_rates.append(float(y[test].mean()))
    return random_rates, stratified_rates, empty


def spread(values) -> dict:
    """Mean, standard deviation, minimum and maximum, rounded for reporting."""
    values = np.asarray(values, dtype=float)
    return {
        "mean": round(float(values.mean()), 4),
        "sd": round(float(values.std()), 4),
        "min": round(float(values.min()), 4),
        "max": round(float(values.max()), 4),
    }


# --------------------------------------------------------------------------
# 3. Groups: when the row is not the unit
# --------------------------------------------------------------------------


def grouped_dataset(n_people: int = 50, rows_each: int = 20, seed: int = 5):
    """Twenty rows per person, and the label is a property of the PERSON.

    There is nothing generalisable here at all: each person's label is a
    coin flip, so a model can only score above chance on a person it has
    already seen. Which is exactly what a row-wise split hands it.
    """
    rng = np.random.default_rng(seed)
    groups = np.repeat(np.arange(n_people), rows_each)
    person_signature = rng.normal(size=(n_people, 4)) * 2.0
    X = person_signature[groups] + rng.normal(size=(n_people * rows_each, 4)) * 0.3
    person_label = rng.integers(0, 2, size=n_people)
    return X, person_label[groups], groups


def rowwise_vs_group_split(X, y, groups, splits: int = 20, test_size: float = 0.25):
    """Score a 1-NN under a row-wise split and under a group-aware one."""
    rowwise = []
    grouped = []
    for seed in range(splits):
        train, test = train_test_split(np.arange(len(y)), test_size=test_size, random_state=seed)
        model = KNeighborsClassifier(1).fit(X[train], y[train])
        rowwise.append(accuracy(y[test], model.predict(X[test])))

        splitter = GroupShuffleSplit(n_splits=1, test_size=test_size, random_state=seed)
        train_g, test_g = next(splitter.split(X, y, groups))
        model_g = KNeighborsClassifier(1).fit(X[train_g], y[train_g])
        grouped.append(accuracy(y[test_g], model_g.predict(X[test_g])))
    return float(np.mean(rowwise)), float(np.mean(grouped))


def groups_shared_between_halves(groups, test_size: float = 0.25, seed: int = 0) -> int:
    """How many people appear in BOTH halves of a row-wise split."""
    train, test = train_test_split(np.arange(len(groups)), test_size=test_size, random_state=seed)
    return len(set(groups[train].tolist()) & set(groups[test].tolist()))


# --------------------------------------------------------------------------
# 4. Time: when the data has a direction
# --------------------------------------------------------------------------


def regime_series(length: int = 1200, n_regimes: int = 6, seed: int = 9):
    """A series in which the rule mapping features to labels changes.

    Six regimes, each with its own randomly drawn linear rule. A random
    shuffle scatters every regime across both halves, so the model always
    has neighbours from the same regime. A chronological split asks it to
    predict a regime it has never seen -- which is what deployment does.
    """
    rng = np.random.default_rng(seed)
    X = rng.normal(size=(length, 4))
    y = np.empty(length, dtype=int)
    block = length // n_regimes
    for r in range(n_regimes):
        rule = rng.normal(size=4)
        window = slice(r * block, (r + 1) * block)
        y[window] = (X[window] @ rule > 0).astype(int)
    return X, y


def shuffled_vs_chronological(X, y, splits: int = 10, test_size: float = 0.25):
    """Score a 5-NN under a shuffled split and under a chronological one.

    Returns ``(mean_shuffled, chronological, majority_baseline_on_the_tail)``.
    """
    shuffled = []
    for seed in range(splits):
        train, test = train_test_split(
            np.arange(len(y)), test_size=test_size, random_state=seed, shuffle=True
        )
        model = KNeighborsClassifier(5).fit(X[train], y[train])
        shuffled.append(accuracy(y[test], model.predict(X[test])))

    cut = int(len(y) * (1.0 - test_size))
    model = KNeighborsClassifier(5).fit(X[:cut], y[:cut])
    chronological = accuracy(y[cut:], model.predict(X[cut:]))
    tail = y[cut:]
    baseline = float(max(tail.mean(), 1.0 - tail.mean()))
    return float(np.mean(shuffled)), chronological, baseline


def temporal_inflation_over_constructions(constructions: int = 20, splits: int = 10):
    """The same comparison over many independently generated series.

    One construction is an anecdote. Reporting the seed that gave the
    largest gap would be the forking-paths problem in a lesson about not
    doing that -- so this returns the whole distribution.
    """
    rows = []
    for seed in range(constructions):
        X, y = regime_series(seed=seed)
        shuffled, chronological, baseline = shuffled_vs_chronological(X, y, splits=splits)
        rows.append((seed, shuffled, chronological, baseline, shuffled - chronological))
    return rows


# --------------------------------------------------------------------------
# 5. One holdout, or many folds
# --------------------------------------------------------------------------


def weak_signal_dataset(n: int = 400, seed: int = 21):
    """A real but modest relationship, so the estimate has something to vary around."""
    rng = np.random.default_rng(seed)
    X = rng.normal(size=(n, 4))
    y = (X[:, 0] + X[:, 1] * 0.6 + rng.normal(size=n) * 1.2 > 0).astype(int)
    return X, y


def holdout_vs_cross_validation(X, y, repeats: int = 200, test_size: float = 0.25, folds: int = 5):
    """The spread of a single-holdout estimate against a k-fold estimate.

    Same data, same model. The only thing that changes between repeats is
    which rows landed where.
    """
    holdout = []
    cross = []
    for seed in range(repeats):
        train, test = train_test_split(
            np.arange(len(y)), test_size=test_size, random_state=seed, stratify=y
        )
        model = LogisticRegression(max_iter=1000).fit(X[train], y[train])
        holdout.append(accuracy(y[test], model.predict(X[test])))

        cross.append(
            float(
                np.mean(
                    cross_val_score(
                        LogisticRegression(max_iter=1000),
                        X,
                        y,
                        cv=StratifiedKFold(folds, shuffle=True, random_state=seed),
                    )
                )
            )
        )
    return holdout, cross


# --------------------------------------------------------------------------
# 6. How big does the test set need to be?
# --------------------------------------------------------------------------


def test_size_table(sizes, p: float = 0.85, draws: int = 20000, seed: int = 33):
    """Predicted and measured standard error of an accuracy, by test-set size.

    Rows are ``(n, theoretical_se, measured_sd, half_width_of_95_interval)``.
    The theory is Day 117's; this is it arriving where the decisions are.
    """
    rows = []
    for n in sizes:
        theory = proportion_standard_error(p, n)
        sample = np.random.default_rng(seed).binomial(n, p, size=draws) / n
        rows.append((n, round(theory, 4), round(float(sample.std()), 4), round(1.96 * theory, 4)))
    return rows


def rows_needed_for_precision(p: float, half_width: float) -> int:
    """Smallest test set whose 95 percent interval is no wider than requested."""
    n = int(np.ceil(p * (1.0 - p) * (1.96 / half_width) ** 2))
    return n


# --------------------------------------------------------------------------
# 7. The rule, made checkable
# --------------------------------------------------------------------------


class TestSetTouchedTwice(RuntimeError):
    """Raised when the test set is evaluated against more than once."""


class GatedTestSet:
    """A test set that permits exactly one evaluation, then refuses.

    Not a substitute for discipline. It is the discipline made mechanical,
    in the same spirit as Day 143's stage contract: a rule that lives in
    code is a rule somebody can check.
    """

    def __init__(self, X, y):
        self._X = X
        self._y = y
        self.evaluations = 0

    def evaluate(self, model) -> float:
        if self.evaluations >= 1:
            raise TestSetTouchedTwice(
                "the test set has already been used once; any further score is a "
                "validation score, not a test score"
            )
        self.evaluations += 1
        return accuracy(self._y, model.predict(self._X))

"""The reference solutions: what each way of splitting a dataset costs.

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

from sklearn.linear_model import LogisticRegression

import splits_lib as s


@pytest.fixture(scope="module")
def rare():
    return s.rare_class_dataset()


@pytest.fixture(scope="module")
def grouped():
    return s.grouped_dataset()


@pytest.fixture(scope="module")
def weak():
    return s.weak_signal_dataset()


# --- 1. Why three sets and not two --------------------------------------


def test_01_picking_the_best_of_k_inflates_the_score_you_picked_it_by():
    rows = s.selection_bias_curve([1, 2, 5, 10, 25, 50, 100, 500, 1000])
    assert rows == [
        (1, 0.4984, 0.5011, -0.0028),
        (2, 0.5115, 0.4992, 0.0123),
        (5, 0.5256, 0.5005, 0.0251),
        (10, 0.5331, 0.4999, 0.0332),
        (25, 0.5436, 0.5009, 0.0427),
        (50, 0.5508, 0.5014, 0.0493),
        (100, 0.5567, 0.4992, 0.0575),
        (500, 0.5682, 0.4978, 0.0704),
        (1000, 0.572, 0.4992, 0.0728),
    ]
    validation = [v for _k, v, _t, _o in rows]
    # The validation score climbs with every extra candidate considered.
    assert all(a < b for a, b in zip(validation, validation[1:]))
    # Every candidate has exactly zero skill: these are coin flips.
    assert rows[-1][1] > rows[0][1] + 0.07


def test_01b_the_test_set_is_the_control_and_stays_at_chance():
    rows = s.selection_bias_curve([1, 10, 100, 1000])
    test_scores = [t for _k, _v, t, _o in rows]
    # Never selected on, therefore never inflated -- at any K.
    for score in test_scores:
        assert abs(score - 0.5) < 0.005
    assert max(test_scores) - min(test_scores) < 0.005
    # This is what makes the validation column mean something.
    validation = [v for _k, v, _t, _o in rows]
    assert validation[-1] - validation[0] > 0.07


def test_01c_the_optimism_is_the_expected_maximum_of_k_noise_draws():
    """And the usual closed-form approximation overestimates it."""
    standard_error = s.proportion_standard_error(0.5, 500)
    assert round(standard_error, 4) == 0.0224

    rows = s.selection_bias_curve([2, 5, 10, 25, 50, 100, 500, 1000])
    for k, _validation, _test, optimism in rows:
        in_errors = optimism / standard_error
        simulated = s.expected_max_of_normals(k)
        approximation = s.sqrt_two_log_k(k)
        # Measurement tracks the simulated expected maximum closely.
        assert abs(in_errors - simulated) < 0.2, (k, in_errors, simulated)
        # The sqrt(2 ln K) asymptotic is above it at every K tried here.
        assert approximation > simulated

    # Concretely, at K = 100: 2.57 standard errors measured, 2.50 expected,
    # against the approximation's 3.03.
    optimism_100 = dict((k, o) for k, _v, _t, o in rows)[100]
    assert round(optimism_100 / standard_error, 2) == 2.57
    assert round(s.expected_max_of_normals(100), 2) == 2.5
    assert round(s.sqrt_two_log_k(100), 2) == 3.03


# --- 2. Stratification ---------------------------------------------------


def test_02_a_random_split_of_a_rare_class_sometimes_has_no_positives(rare):
    X, y = rare
    assert round(float(y.mean()), 4) == 0.05
    random_rates, stratified_rates, empty = s.split_positive_rates(X, y)
    assert s.spread(random_rates) == {"mean": 0.0504, "sd": 0.0265, "min": 0.0, "max": 0.16}
    assert s.spread(stratified_rates) == {"mean": 0.05, "sd": 0.01, "min": 0.04, "max": 0.06}
    # Twenty-one splits in five hundred produced a test set with no
    # positives at all, on which recall is undefined.
    assert empty == 21
    assert 0 in [round(r, 4) for r in random_rates]
    assert 0 not in [round(r, 4) for r in stratified_rates]


def test_02b_stratifying_shrinks_the_spread_without_changing_the_mean(rare):
    X, y = rare
    random_rates, stratified_rates, _empty = s.split_positive_rates(X, y)
    random = s.spread(random_rates)
    stratified = s.spread(stratified_rates)
    assert abs(random["mean"] - stratified["mean"]) < 0.001
    assert stratified["sd"] < random["sd"]
    assert round(random["sd"] / stratified["sd"], 2) == 2.65
    # The random split's worst case is three times the population rate.
    assert round(random["max"] / float(y.mean()), 2) == 3.2


# --- 3. Groups -----------------------------------------------------------


def test_03_splitting_rows_when_the_unit_is_a_person_invents_fifty_six_points(grouped):
    X, y, groups = grouped
    rowwise, group_aware = s.rowwise_vs_group_split(X, y, groups)
    assert round(rowwise, 4) == 0.976
    assert round(group_aware, 4) == 0.4112
    assert round(rowwise - group_aware, 4) == 0.5648
    # There is nothing generalisable in this data at all: each person's
    # label is a coin flip, so the group-aware score is chance, and low.
    assert group_aware < 0.5 < rowwise


def test_03b_every_single_person_appears_in_both_halves(grouped):
    _X, _y, groups = grouped
    shared = s.groups_shared_between_halves(groups)
    assert shared == 50
    assert len(set(groups.tolist())) == 50
    # Twenty rows each: a random quarter cannot miss anybody.
    assert len(groups) == 1000


# --- 4. Time -------------------------------------------------------------


def test_04_a_shuffled_split_beats_a_chronological_one_every_single_time():
    rows = s.temporal_inflation_over_constructions()
    assert len(rows) == 20
    inflation = [gap for _seed, _sh, _ch, _base, gap in rows]
    # The direction is universal across twenty independent constructions.
    assert all(gap > 0 for gap in inflation)
    assert sum(1 for gap in inflation if gap > 0) == 20


def test_04b_but_the_size_of_the_effect_varies_by_an_order_of_magnitude():
    rows = s.temporal_inflation_over_constructions()
    shuffled = float(np.mean([sh for _s, sh, _c, _b, _g in rows]))
    chronological = float(np.mean([ch for _s, _sh, ch, _b, _g in rows]))
    baseline = float(np.mean([b for _s, _sh, _ch, b, _g in rows]))
    inflation = [gap for _s, _sh, _ch, _b, gap in rows]
    assert round(shuffled, 4) == 0.5961
    assert round(chronological, 4) == 0.5233
    assert round(baseline, 4) == 0.5235
    assert round(float(np.mean(inflation)), 4) == 0.0728
    assert round(float(np.std(inflation)), 4) == 0.0596
    assert round(min(inflation), 4) == 0.016
    assert round(max(inflation), 4) == 0.2557
    # Sixteen times between the smallest and largest effect: reporting the
    # largest would be the forking-paths problem in a lesson against it.
    assert round(max(inflation) / min(inflation), 1) == 16.0
    # The honest verdict: chronologically, the model has learned nothing.
    assert abs(chronological - baseline) < 0.005
    assert shuffled > baseline


# --- 5. One holdout, or many folds --------------------------------------


def test_05_one_holdout_swings_nineteen_points_on_identical_data(weak):
    X, y = weak
    holdout, cross = s.holdout_vs_cross_validation(X, y)
    assert s.spread(holdout) == {"mean": 0.7519, "sd": 0.0381, "min": 0.66, "max": 0.85}
    assert s.spread(cross) == {"mean": 0.7546, "sd": 0.0061, "min": 0.7375, "max": 0.77}
    # Same data, same model. Only which rows landed where changed.
    assert round(max(holdout) - min(holdout), 4) == 0.19
    assert round(max(cross) - min(cross), 4) == 0.0325


def test_05b_cross_validation_is_six_times_steadier_for_the_same_data(weak):
    X, y = weak
    holdout, cross = s.holdout_vs_cross_validation(X, y)
    ratio = float(np.std(holdout) / np.std(cross))
    assert round(ratio, 4) == 6.2344
    # It estimates the same thing -- the means agree to within 0.003.
    assert abs(float(np.mean(holdout)) - float(np.mean(cross))) < 0.003


# --- 6. How big must the test set be? -----------------------------------


def test_06_the_standard_error_formula_predicts_the_measured_spread():
    rows = s.test_size_table([50, 100, 200, 500, 1000, 5000])
    assert rows == [
        (50, 0.0505, 0.0505, 0.099),
        (100, 0.0357, 0.0357, 0.07),
        (200, 0.0252, 0.0254, 0.0495),
        (500, 0.016, 0.016, 0.0313),
        (1000, 0.0113, 0.0112, 0.0221),
        (5000, 0.005, 0.0051, 0.0099),
    ]
    for _n, theory, measured, _half in rows:
        assert abs(theory - measured) <= 0.0002
    # Four times the rows halves the error, not quarters it.
    by_n = {n: theory for n, theory, _m, _h in rows}
    assert round(by_n[200] / by_n[50], 2) == 0.5
    assert round(by_n[5000] / by_n[500], 2) == 0.31


def test_06b_a_hundred_row_test_set_cannot_resolve_a_five_point_difference():
    half_width = s.test_size_table([100])[0][3]
    assert half_width == 0.07
    # Plus or minus seven points: two models five points apart are
    # indistinguishable on it.
    assert half_width > 0.05
    assert s.rows_needed_for_precision(0.85, 0.02) == 1225
    assert s.rows_needed_for_precision(0.85, 0.01) == 4899


# --- 7. The rule, made mechanical ---------------------------------------


def test_07_the_test_set_permits_exactly_one_evaluation(weak):
    X, y = weak
    model = LogisticRegression(max_iter=1000).fit(X, y)
    gate = s.GatedTestSet(X, y)
    assert gate.evaluations == 0
    first = gate.evaluate(model)
    assert round(first, 4) == 0.7575
    assert gate.evaluations == 1
    with pytest.raises(s.TestSetTouchedTwice) as excinfo:
        gate.evaluate(model)
    assert "validation score" in str(excinfo.value)
    # And the counter does not advance on a refused attempt.
    assert gate.evaluations == 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 splits_lib as s


def test_the_standard_error_formula_behaves_as_it_should():
    # Maximal at p = 0.5, and shrinking like one over root n.
    assert s.proportion_standard_error(0.5, 100) > s.proportion_standard_error(0.85, 100)
    assert s.proportion_standard_error(0.5, 100) > s.proportion_standard_error(0.5, 400)
    quartered = s.proportion_standard_error(0.5, 400) / s.proportion_standard_error(0.5, 100)
    assert round(quartered, 4) == 0.5
    # A certain outcome has no sampling error at all.
    assert s.proportion_standard_error(1.0, 100) == 0.0


def test_the_grouped_dataset_really_is_grouped():
    X, y, groups = s.grouped_dataset(n_people=10, rows_each=5, seed=1)
    assert X.shape == (50, 4) and y.shape == (50,) and groups.shape == (50,)
    assert len(set(groups.tolist())) == 10
    # Every row belonging to one person carries that person's single label.
    for person in range(10):
        member_labels = set(y[groups == person].tolist())
        assert len(member_labels) == 1


def test_the_regime_series_really_changes_its_rule():
    X, y = s.regime_series(length=600, n_regimes=3, seed=2)
    assert X.shape == (600, 4) and y.shape == (600,)
    # Each regime is internally consistent but the regimes disagree: a
    # linear model fitted on one block should do worse on another.
    from sklearn.linear_model import LogisticRegression

    first = LogisticRegression(max_iter=1000).fit(X[:200], y[:200])
    own = s.accuracy(y[:200], first.predict(X[:200]))
    other = s.accuracy(y[400:], first.predict(X[400:]))
    assert own > 0.9
    assert other < own - 0.2


def test_the_gated_test_set_counts_and_refuses():
    class AlwaysZero:
        def predict(self, X):
            return np.zeros(len(X), dtype=int)

    y = np.array([0, 0, 0, 1])
    gate = s.GatedTestSet(np.zeros((4, 2)), y)
    assert gate.evaluate(AlwaysZero()) == 0.75
    with pytest.raises(s.TestSetTouchedTwice):
        gate.evaluate(AlwaysZero())
    # A fresh gate is a fresh budget; the class holds no global state.
    assert s.GatedTestSet(np.zeros((4, 2)), y).evaluate(AlwaysZero()) == 0.75

lesson_id: D144
day: 144
kind: guided-build
languages:
  - python
  - bash
setup_commands:
  - cd labs/sections/machine-learning/day-144-train-validation-and-test-splits
  - 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: 65
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_splits_lib.py are solved in both directories; the fourteen exercise stubs in
  starter/test_splits_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 round(rowwise -
  group_aware, 4) == 0.5648` to 0.0, confirms a non-zero exit naming the failing test,
  and removes the scratch directory. Separately, by hand, `assert empty == 21` was
  changed to 999 in examples/test_splits_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. (1) WHY THREE SETS: K candidates that are literal
  coin flips, each scored on a 500-row validation set and a 500-row test set, averaged
  over 400 replications. Picking the best on validation gives 0.4984, 0.5115, 0.5256,
  0.5331, 0.5436, 0.5508, 0.5567, 0.5682 and 0.5720 at K = 1, 2, 5, 10, 25, 50, 100, 500
  and 1000 -- climbing steadily on candidates with exactly zero skill -- while that same
  candidate's test score never moves from chance, staying within 0.003 of 0.5 at every
  K. The test column is the control and it is what makes the validation column mean
  anything. (2) The optimism is not a vague warning but a computable quantity: expressed
  in standard errors (0.0224 at n=500) it reads 0.55, 1.12, 1.48, 1.91, 2.20, 2.57, 3.15
  and 3.26, tracking the simulated expected maximum of K standard normals (0.55, 1.16,
  1.54, 1.97, 2.25, 2.50, 3.03, 3.24) to within 0.2 at every K -- while the familiar
  closed form sqrt(2 ln K) sits ABOVE the truth at every K here (1.18, 1.79, 2.15, 2.54,
  2.80, 3.03, 3.53, 3.72). (3) Stratification: on a 5-percent positive class, 500 random
  splits give test positive rates with sd 0.0265 against 0.0100 stratified, and 21 of
  those 500 random splits produced a test half containing NO positives at all, on which
  recall is undefined. (4) Groups: 50 people, 20 rows each, each person's label a coin
  flip so nothing is generalisable. A row-wise random split scores 0.9760; a group-aware
  split scores 0.4112; 0.5648 of accuracy invented. All 50 people appear in both halves
  of the row-wise split -- with 20 rows each, a random quarter cannot miss anybody. (5)
  Time: across 20 independently generated six-regime series, a shuffled split beat a
  chronological one in 20 of 20, with shuffled mean 0.5961 against chronological 0.5233
  and a majority baseline of 0.5235 -- so chronologically the model has learned nothing.
  (6) One holdout against 5-fold on identical data and an identical model, over 200
  seeds: the holdout has sd 0.0381 and spans 0.66 to 0.85, a swing of 0.19, while 5-fold
  has sd 0.0061 and spans 0.7375 to 0.77; cross-validation is 6.2344 times steadier and
  the two means agree to within 0.003. (7) Test-set size: the predicted standard error
  sqrt(p(1-p)/n) matches 20000 measured draws to within 0.0002 at n = 50, 100, 200, 500,
  1000 and 5000, so a 100-row test set carries a 95-percent half-width of 0.07 and
  cannot resolve a five-point difference; 1225 rows are needed for plus or minus 0.02 at
  an accuracy of 0.85, and 4899 for plus or minus 0.01. (8) A GatedTestSet permits
  exactly one evaluation, returns 0.7575, and raises TestSetTouchedTwice on the second
  with a message naming what the second number really is -- without advancing its
  counter on the refused attempt. THREE HONESTY CALLS. FIRST, and the most important:
  the temporal effect's SIZE varies by a factor of sixteen across constructions, from
  +0.0160 to +0.2557. The first construction tried gave +0.1428 and quoting it alone
  would have been the forking-paths problem inside a lab against exactly that, so
  exercise 4 asserts only the direction, which held 20 times out of 20, and exercise 4b
  asserts the whole distribution -- mean 0.0728, sd 0.0596, min and max -- rather than a
  single headline figure. Not every split mistake costs the same: group leakage cost 56
  points and temporal leakage cost 7 on average. SECOND: the group-aware score of 0.4112
  is BELOW chance and this is not anti-learning -- about twelve people land in that test
  half, each contributing twenty identical labels, so the estimate is roughly a dozen
  coin flips and it wanders; the lab asserts the structural claim (below 0.5 while
  row-wise is far above it) rather than the value. THIRD: the textbook sqrt(2 ln K)
  approximation is reported as overestimating the measured optimism at every K tried,
  because that is what was measured; the lab asserts the inequality rather than treating
  the closed form as the truth. Harness check 8 re-runs group leakage at five dataset
  seeds, selection optimism at a different replication count, and stratification at
  three seeds, so every direction the lesson claims is confirmed beyond the seed it
  quotes.

# 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

Almost every number in this lab is an average over seeded draws from
`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 sampled figure would move.

What does not depend on the pins: the standard-error formula in exercise
6, which is arithmetic; the direction of every result — group leakage
inflates, stratification narrows, shuffling beats chronology, selecting on
a set inflates that set's score; and the structural facts, such as all
fifty people appearing in both halves of a row-wise split. Harness check 8
re-runs three of those directions at 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 144 lab brief — Three Sets, and Why

Everybody knows you hold out a test set. Rather fewer people can say why
there are supposed to be *three* sets rather than two, and almost nobody
has seen the number that justifies the third one.

This lab measures it.

## The claim you are here to measure

> A validation set you select on is a set you have fitted to.

Exercise 1 does it with candidates that have exactly zero skill. Each one
is a coin flip: a fixed vector of random predictions, scored on a 500-row
validation set and a 500-row test set. Pick whichever scores best on
validation, then look at what it scores on test:

| candidates considered | best validation | its test score | optimism |
| --- | --- | --- | --- |
| 1 | 0.4984 | 0.5011 | −0.0028 |
| 10 | 0.5331 | 0.4999 | +0.0332 |
| 100 | 0.5567 | 0.4992 | +0.0575 |
| 1000 | 0.5720 | 0.4992 | +0.0728 |

Read the test column first. **It never moves.** It sits at chance for
every K, because it was never selected on. That is the control, and it is
what makes the validation column mean something.

Now read the validation column. It climbs, steadily, all the way to
0.5720 — on candidates that are coin flips. Try a thousand things and the
best of them will look seven points better than chance, every time,
whether or not any of them is any good.

That is why there are three sets. Not convention. Arithmetic.

## The part that is genuinely satisfying

Exercise 1c checks the optimism against theory. Express it in standard
errors — the standard error of an accuracy on 500 rows is 0.0224 — and
compare it to the expected maximum of K standard normal draws:

| K | measured, in SEs | E max of K normals | sqrt(2 ln K) |
| --- | --- | --- | --- |
| 10 | 1.48 | 1.54 | 2.15 |
| 100 | 2.57 | 2.50 | 3.03 |
| 1000 | 3.26 | 3.24 | 3.72 |

The measurement tracks the simulated expectation closely — and the
familiar closed-form `sqrt(2 ln K)` sits above the truth at every K here.
It is an asymptotic, and it is loose at any K you will actually use. The
lab asserts the inequality rather than a gap size.

## The four ways a split goes wrong

| # | The mistake | What it costs, measured |
| --- | --- | --- |
| 2 | not stratifying a rare class | 21 of 500 random splits had a test half with **no positives at all** |
| 3 | splitting rows when the unit is a person | **+0.5648** — 0.9760 against 0.4112, and all 50 people were in both halves |
| 4 | shuffling data that has a direction in time | +0.0728 on average, and shuffling won in **20 of 20** constructions |
| 5 | reading a trend off one holdout | one holdout swung **0.19** across seeds; 5-fold swung 0.0325 |

Exercise 3 is the one that should alarm you. Fifty people, twenty rows
each, and each person's label is a coin flip — there is nothing
generalisable in that dataset whatsoever. A row-wise random split reports
**97.6 percent accuracy**. A group-aware split reports 0.4112, which is
chance. Fifty-six points, and the mechanism is one line: with twenty rows
each, a random quarter cannot miss anybody, so every test person is
already in training.

## The honesty call in exercise 4

The temporal effect is real: shuffling beat chronology in all twenty
constructions. But its **size varies by a factor of sixteen**, from +0.016
to +0.2557.

The first construction tried while building this lab gave +0.1428.
Reporting that one number would have been the forking-paths problem inside
a lab about not committing it. So exercise 4 splits in two: 4 asserts the
direction, which held every time, and 4b asserts the whole distribution —
mean, standard deviation, minimum and maximum — and the fact that the
chronological score is statistically indistinguishable from the majority
baseline.

Not every split mistake costs the same. Group leakage cost fifty-six
points. Temporal leakage, here, cost seven on average. Both are real; only
one is an emergency.

## 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_splits_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 `splits_lib.py`, `test_splits_lib.py` and `test_splits_claims.py`;
pytest aborts on the module-name collision. Run them separately, always.

## And the rule, made mechanical

Exercise 7 wraps the test set in a `GatedTestSet` that permits exactly one
evaluation and raises `TestSetTouchedTwice` on the second, with a message
saying what the second number actually is: a validation score.

That is not a substitute for discipline. It is the discipline made
mechanical, in the same spirit as Day 143's stage contract — a rule that
lives in code is a rule somebody can check.

"""Splitting, measured: what each way of cutting a dataset actually buys.

Three sets, not two, and the reason is arithmetic rather than convention.
A validation set you select on is a set you have fitted to, and this
module measures the resulting optimism directly -- it turns out to be
exactly the expected maximum of K noise draws, which is a quantity you can
compute.

The rest measures the four ways a split goes wrong: not stratifying when
the class is rare, splitting rows when the unit is a person, splitting
randomly when the data has a direction in time, and reading a trend off
one holdout when a holdout's own spread is wider than the trend.

Everything here is deterministic given a seed.
"""

from __future__ import annotations

import numpy as np

from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import (
    GroupShuffleSplit,
    StratifiedKFold,
    StratifiedShuffleSplit,
    cross_val_score,
    train_test_split,
)
from sklearn.neighbors import KNeighborsClassifier


def accuracy(y_true, y_pred) -> float:
    return float(np.mean(np.asarray(y_true) == np.asarray(y_pred)))


# --------------------------------------------------------------------------
# 1. Why three sets: selecting on a set is fitting to it
# --------------------------------------------------------------------------


def selection_replicate(k_candidates: int, n: int = 500, seed: int = 0):
    """One replication of "pick the best of K, then check it on a fresh set".

    Each candidate is a fixed prediction vector over 2n rows, split into an
    n-row validation set and an n-row test set. Every candidate has exactly
    zero skill by construction, so any validation score above 0.5 is noise
    and any *selected* validation score above 0.5 is noise you chose.

    Returns ``(best_validation_score, that_candidate_s_test_score)``.
    """
    rng = np.random.default_rng(seed)
    y = rng.integers(0, 2, size=2 * n)
    predictions = rng.integers(0, 2, size=(k_candidates, 2 * n))
    correct = predictions == y
    validation = correct[:, :n].mean(axis=1)
    test = correct[:, n:].mean(axis=1)
    winner = int(np.argmax(validation))
    return float(validation[winner]), float(test[winner])


def selection_bias_curve(k_values, replications: int = 400, n: int = 500):
    """Mean selected-validation score and its test score, for each K.

    Returns rows of ``(k, mean_validation, mean_test, optimism)``. The test
    column is the control: it must stay at chance for every K, because the
    test set was never selected on.
    """
    rows = []
    for k in k_values:
        pairs = [selection_replicate(k, n=n, seed=r) for r in range(replications)]
        validation = float(np.mean([v for v, _t in pairs]))
        test = float(np.mean([t for _v, t in pairs]))
        rows.append((k, round(validation, 4), round(test, 4), round(validation - test, 4)))
    return rows


def proportion_standard_error(p: float, n: int) -> float:
    """The standard error of an accuracy estimated on n rows."""
    return float(np.sqrt(p * (1.0 - p) / n))


def expected_max_of_normals(k: int, draws: int = 20000, seed: int = 7) -> float:
    """E of the maximum of k standard normals, by simulation.

    This is the quantity the selection optimism should equal, once the
    optimism is expressed in standard errors. Simulated rather than
    approximated, because the usual closed form overestimates it -- which
    exercise 1c measures.
    """
    rng = np.random.default_rng(seed)
    return float(np.mean(np.max(rng.standard_normal((draws, k)), axis=1)))


def sqrt_two_log_k(k: int) -> float:
    """The textbook asymptotic for the expected maximum of k normals."""
    return 0.0 if k <= 1 else float(np.sqrt(2.0 * np.log(k)))


# --------------------------------------------------------------------------
# 2. Stratification: a rare class and a small test set
# --------------------------------------------------------------------------


def rare_class_dataset(n: int = 200, rate: float = 0.05, seed: int = 144):
    """A dataset with a genuinely rare positive class."""
    rng = np.random.default_rng(seed)
    y = np.zeros(n, dtype=int)
    y[: int(n * rate)] = 1
    X = rng.normal(size=(n, 3))
    return X, y


def split_positive_rates(X, y, splits: int = 500, test_size: float = 0.25):
    """Positive rate in the test half, over many random and stratified splits.

    Returns ``(random_rates, stratified_rates, random_splits_with_no_positives)``.
    """
    random_rates = []
    stratified_rates = []
    empty = 0
    for seed in range(splits):
        _x_tr, _x_te, _y_tr, y_te = train_test_split(X, y, test_size=test_size, random_state=seed)
        random_rates.append(float(y_te.mean()))
        if y_te.sum() == 0:
            empty += 1
        splitter = StratifiedShuffleSplit(n_splits=1, test_size=test_size, random_state=seed)
        _train, test = next(splitter.split(X, y))
        stratified_rates.append(float(y[test].mean()))
    return random_rates, stratified_rates, empty


def spread(values) -> dict:
    """Mean, standard deviation, minimum and maximum, rounded for reporting."""
    values = np.asarray(values, dtype=float)
    return {
        "mean": round(float(values.mean()), 4),
        "sd": round(float(values.std()), 4),
        "min": round(float(values.min()), 4),
        "max": round(float(values.max()), 4),
    }


# --------------------------------------------------------------------------
# 3. Groups: when the row is not the unit
# --------------------------------------------------------------------------


def grouped_dataset(n_people: int = 50, rows_each: int = 20, seed: int = 5):
    """Twenty rows per person, and the label is a property of the PERSON.

    There is nothing generalisable here at all: each person's label is a
    coin flip, so a model can only score above chance on a person it has
    already seen. Which is exactly what a row-wise split hands it.
    """
    rng = np.random.default_rng(seed)
    groups = np.repeat(np.arange(n_people), rows_each)
    person_signature = rng.normal(size=(n_people, 4)) * 2.0
    X = person_signature[groups] + rng.normal(size=(n_people * rows_each, 4)) * 0.3
    person_label = rng.integers(0, 2, size=n_people)
    return X, person_label[groups], groups


def rowwise_vs_group_split(X, y, groups, splits: int = 20, test_size: float = 0.25):
    """Score a 1-NN under a row-wise split and under a group-aware one."""
    rowwise = []
    grouped = []
    for seed in range(splits):
        train, test = train_test_split(np.arange(len(y)), test_size=test_size, random_state=seed)
        model = KNeighborsClassifier(1).fit(X[train], y[train])
        rowwise.append(accuracy(y[test], model.predict(X[test])))

        splitter = GroupShuffleSplit(n_splits=1, test_size=test_size, random_state=seed)
        train_g, test_g = next(splitter.split(X, y, groups))
        model_g = KNeighborsClassifier(1).fit(X[train_g], y[train_g])
        grouped.append(accuracy(y[test_g], model_g.predict(X[test_g])))
    return float(np.mean(rowwise)), float(np.mean(grouped))


def groups_shared_between_halves(groups, test_size: float = 0.25, seed: int = 0) -> int:
    """How many people appear in BOTH halves of a row-wise split."""
    train, test = train_test_split(np.arange(len(groups)), test_size=test_size, random_state=seed)
    return len(set(groups[train].tolist()) & set(groups[test].tolist()))


# --------------------------------------------------------------------------
# 4. Time: when the data has a direction
# --------------------------------------------------------------------------


def regime_series(length: int = 1200, n_regimes: int = 6, seed: int = 9):
    """A series in which the rule mapping features to labels changes.

    Six regimes, each with its own randomly drawn linear rule. A random
    shuffle scatters every regime across both halves, so the model always
    has neighbours from the same regime. A chronological split asks it to
    predict a regime it has never seen -- which is what deployment does.
    """
    rng = np.random.default_rng(seed)
    X = rng.normal(size=(length, 4))
    y = np.empty(length, dtype=int)
    block = length // n_regimes
    for r in range(n_regimes):
        rule = rng.normal(size=4)
        window = slice(r * block, (r + 1) * block)
        y[window] = (X[window] @ rule > 0).astype(int)
    return X, y


def shuffled_vs_chronological(X, y, splits: int = 10, test_size: float = 0.25):
    """Score a 5-NN under a shuffled split and under a chronological one.

    Returns ``(mean_shuffled, chronological, majority_baseline_on_the_tail)``.
    """
    shuffled = []
    for seed in range(splits):
        train, test = train_test_split(
            np.arange(len(y)), test_size=test_size, random_state=seed, shuffle=True
        )
        model = KNeighborsClassifier(5).fit(X[train], y[train])
        shuffled.append(accuracy(y[test], model.predict(X[test])))

    cut = int(len(y) * (1.0 - test_size))
    model = KNeighborsClassifier(5).fit(X[:cut], y[:cut])
    chronological = accuracy(y[cut:], model.predict(X[cut:]))
    tail = y[cut:]
    baseline = float(max(tail.mean(), 1.0 - tail.mean()))
    return float(np.mean(shuffled)), chronological, baseline


def temporal_inflation_over_constructions(constructions: int = 20, splits: int = 10):
    """The same comparison over many independently generated series.

    One construction is an anecdote. Reporting the seed that gave the
    largest gap would be the forking-paths problem in a lesson about not
    doing that -- so this returns the whole distribution.
    """
    rows = []
    for seed in range(constructions):
        X, y = regime_series(seed=seed)
        shuffled, chronological, baseline = shuffled_vs_chronological(X, y, splits=splits)
        rows.append((seed, shuffled, chronological, baseline, shuffled - chronological))
    return rows


# --------------------------------------------------------------------------
# 5. One holdout, or many folds
# --------------------------------------------------------------------------


def weak_signal_dataset(n: int = 400, seed: int = 21):
    """A real but modest relationship, so the estimate has something to vary around."""
    rng = np.random.default_rng(seed)
    X = rng.normal(size=(n, 4))
    y = (X[:, 0] + X[:, 1] * 0.6 + rng.normal(size=n) * 1.2 > 0).astype(int)
    return X, y


def holdout_vs_cross_validation(X, y, repeats: int = 200, test_size: float = 0.25, folds: int = 5):
    """The spread of a single-holdout estimate against a k-fold estimate.

    Same data, same model. The only thing that changes between repeats is
    which rows landed where.
    """
    holdout = []
    cross = []
    for seed in range(repeats):
        train, test = train_test_split(
            np.arange(len(y)), test_size=test_size, random_state=seed, stratify=y
        )
        model = LogisticRegression(max_iter=1000).fit(X[train], y[train])
        holdout.append(accuracy(y[test], model.predict(X[test])))

        cross.append(
            float(
                np.mean(
                    cross_val_score(
                        LogisticRegression(max_iter=1000),
                        X,
                        y,
                        cv=StratifiedKFold(folds, shuffle=True, random_state=seed),
                    )
                )
            )
        )
    return holdout, cross


# --------------------------------------------------------------------------
# 6. How big does the test set need to be?
# --------------------------------------------------------------------------


def test_size_table(sizes, p: float = 0.85, draws: int = 20000, seed: int = 33):
    """Predicted and measured standard error of an accuracy, by test-set size.

    Rows are ``(n, theoretical_se, measured_sd, half_width_of_95_interval)``.
    The theory is Day 117's; this is it arriving where the decisions are.
    """
    rows = []
    for n in sizes:
        theory = proportion_standard_error(p, n)
        sample = np.random.default_rng(seed).binomial(n, p, size=draws) / n
        rows.append((n, round(theory, 4), round(float(sample.std()), 4), round(1.96 * theory, 4)))
    return rows


def rows_needed_for_precision(p: float, half_width: float) -> int:
    """Smallest test set whose 95 percent interval is no wider than requested."""
    n = int(np.ceil(p * (1.0 - p) * (1.96 / half_width) ** 2))
    return n


# --------------------------------------------------------------------------
# 7. The rule, made checkable
# --------------------------------------------------------------------------


class TestSetTouchedTwice(RuntimeError):
    """Raised when the test set is evaluated against more than once."""


class GatedTestSet:
    """A test set that permits exactly one evaluation, then refuses.

    Not a substitute for discipline. It is the discipline made mechanical,
    in the same spirit as Day 143's stage contract: a rule that lives in
    code is a rule somebody can check.
    """

    def __init__(self, X, y):
        self._X = X
        self._y = y
        self.evaluations = 0

    def evaluate(self, model) -> float:
        if self.evaluations >= 1:
            raise TestSetTouchedTwice(
                "the test set has already been used once; any further score is a "
                "validation score, not a test score"
            )
        self.evaluations += 1
        return accuracy(self._y, model.predict(self._X))

"""Fourteen exercises in what each way of splitting a dataset costs.

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.
`splits_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

from sklearn.linear_model import LogisticRegression  # noqa: F401

import splits_lib as s  # noqa: F401  (you will need it)


@pytest.fixture(scope="module")
def rare():
    return s.rare_class_dataset()


@pytest.fixture(scope="module")
def grouped():
    return s.grouped_dataset()


@pytest.fixture(scope="module")
def weak():
    return s.weak_signal_dataset()


def test_01_picking_the_best_of_k_inflates_the_score_you_picked_it_by():
    pytest.skip(
        "Assert s.selection_bias_curve([1, 2, 5, 10, 25, 50, 100, 500, 1000]) "
        "equals the nine rows in expected-output/measured-values.txt, from "
        "(1, 0.4984, 0.5011, -0.0028) to (1000, 0.572, 0.4992, 0.0728). Then "
        "assert the validation column is strictly increasing in K. Every "
        "candidate is a coin flip with exactly zero skill, so all of that "
        "climb is noise you selected."
    )


def test_01b_the_test_set_is_the_control_and_stays_at_chance():
    pytest.skip(
        "For K in [1, 10, 100, 1000], assert every test score is within "
        "0.005 of 0.5 and that the whole test column spans less than 0.005, "
        "while the validation column climbs by more than 0.07. The test set "
        "was never selected on, so it was never inflated -- and that is what "
        "makes the validation column mean anything."
    )


def test_01c_the_optimism_is_the_expected_maximum_of_k_noise_draws():
    pytest.skip(
        "Assert s.proportion_standard_error(0.5, 500) rounds to 0.0224. For "
        "each K in [2, 5, 10, 25, 50, 100, 500, 1000], divide the optimism "
        "by that standard error and assert it is within 0.2 of "
        "s.expected_max_of_normals(k), and that s.sqrt_two_log_k(k) is "
        "LARGER than the simulated expectation at every K. Then assert the "
        "concrete case: at K=100 the measured optimism is 2.57 standard "
        "errors, the simulated expectation is 2.50, and the textbook "
        "approximation says 3.03."
    )


def test_02_a_random_split_of_a_rare_class_sometimes_has_no_positives(rare):
    pytest.skip(
        "The population positive rate is 0.05. Assert s.split_positive_rates "
        "gives random spread {'mean': 0.0504, 'sd': 0.0265, 'min': 0.0, "
        "'max': 0.16} and stratified spread {'mean': 0.05, 'sd': 0.01, "
        "'min': 0.04, 'max': 0.06}, and that exactly 21 of the 500 random "
        "splits produced a test half with no positives at all. Recall is "
        "undefined on those."
    )


def test_02b_stratifying_shrinks_the_spread_without_changing_the_mean(rare):
    pytest.skip(
        "Assert the two means agree to within 0.001 while the stratified "
        "standard deviation is smaller, with a ratio of exactly 2.65. Then "
        "assert the random split's worst case is 3.2 times the population "
        "positive rate. Stratifying does not change what you are estimating; "
        "it changes how much the estimate wobbles."
    )


def test_03_splitting_rows_when_the_unit_is_a_person_invents_fifty_six_points(grouped):
    pytest.skip(
        "Assert s.rowwise_vs_group_split gives 0.976 for the row-wise split "
        "and 0.4112 for the group-aware one, a gap of 0.5648. Then assert "
        "the group-aware score is below 0.5 and the row-wise one above it. "
        "Each person's label is a coin flip, so there is nothing "
        "generalisable here at all -- and a row-wise split reports 97.6 "
        "percent."
    )


def test_03b_every_single_person_appears_in_both_halves(grouped):
    pytest.skip(
        "Assert s.groups_shared_between_halves returns 50, that there are 50 "
        "distinct people, and that the dataset has 1000 rows. With twenty "
        "rows each, a random quarter cannot miss anybody -- which is the "
        "mechanism behind exercise 3, stated as a count rather than as an "
        "intuition."
    )


def test_04_a_shuffled_split_beats_a_chronological_one_every_single_time():
    pytest.skip(
        "Call s.temporal_inflation_over_constructions(). Assert there are 20 "
        "rows and that the inflation is positive in every single one. The "
        "direction of this effect is universal; the next exercise is about "
        "its size."
    )


def test_04b_but_the_size_of_the_effect_varies_by_an_order_of_magnitude():
    pytest.skip(
        "Assert the means across 20 constructions are shuffled 0.5961, "
        "chronological 0.5233 and baseline 0.5235, with inflation mean "
        "0.0728, sd 0.0596, min 0.016 and max 0.2557 -- a factor of 16.0 "
        "between smallest and largest. Assert the chronological score is "
        "within 0.005 of the baseline while the shuffled one is above it. "
        "Quoting the 0.2557 seed alone would be the forking-paths problem, "
        "in a lab against it."
    )


def test_05_one_holdout_swings_nineteen_points_on_identical_data(weak):
    pytest.skip(
        "Assert s.holdout_vs_cross_validation gives holdout spread {'mean': "
        "0.7519, 'sd': 0.0381, 'min': 0.66, 'max': 0.85} and 5-fold spread "
        "{'mean': 0.7546, 'sd': 0.0061, 'min': 0.7375, 'max': 0.77}. Assert "
        "the holdout range is exactly 0.19 and the cross-validated range "
        "0.0325. Same data, same model: only which rows landed where changed."
    )


def test_05b_cross_validation_is_six_times_steadier_for_the_same_data(weak):
    pytest.skip(
        "Assert the ratio of the two standard deviations is 6.2344, and that "
        "the two means agree to within 0.003. Cross-validation estimates the "
        "same quantity; it just estimates it with far less noise, because "
        "every row serves as test data exactly once."
    )


def test_06_the_standard_error_formula_predicts_the_measured_spread():
    pytest.skip(
        "Assert s.test_size_table([50, 100, 200, 500, 1000, 5000]) matches "
        "the captured rows, and that theory and measurement differ by at "
        "most 0.0002 at every size. Then assert the scaling: quadrupling the "
        "rows halves the error (200 against 50 gives a ratio of 0.5), and "
        "5000 against 500 gives 0.31. This is Day 117's formula arriving "
        "where the decisions are."
    )


def test_06b_a_hundred_row_test_set_cannot_resolve_a_five_point_difference():
    pytest.skip(
        "Assert the 95 percent half-width at n=100 is 0.07, which is wider "
        "than 0.05 -- so two models five points apart are indistinguishable "
        "on it. Then assert s.rows_needed_for_precision(0.85, 0.02) is 1225 "
        "and (0.85, 0.01) is 4899. Decide how big your test set must be "
        "BEFORE you split, from the difference you need to detect."
    )


def test_07_the_test_set_permits_exactly_one_evaluation(weak):
    pytest.skip(
        "Fit LogisticRegression(max_iter=1000) on the weak-signal data, wrap "
        "the data in s.GatedTestSet, and assert the first evaluation is "
        "0.7575 and the counter becomes 1. Then assert a second evaluation "
        "raises s.TestSetTouchedTwice with a message mentioning 'validation "
        "score', and that the counter did NOT advance on the refused "
        "attempt. The gate is not a substitute for discipline; it is the "
        "discipline made mechanical."
    )

"""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 splits_lib as s


def test_the_standard_error_formula_behaves_as_it_should():
    # Maximal at p = 0.5, and shrinking like one over root n.
    assert s.proportion_standard_error(0.5, 100) > s.proportion_standard_error(0.85, 100)
    assert s.proportion_standard_error(0.5, 100) > s.proportion_standard_error(0.5, 400)
    quartered = s.proportion_standard_error(0.5, 400) / s.proportion_standard_error(0.5, 100)
    assert round(quartered, 4) == 0.5
    # A certain outcome has no sampling error at all.
    assert s.proportion_standard_error(1.0, 100) == 0.0


def test_the_grouped_dataset_really_is_grouped():
    X, y, groups = s.grouped_dataset(n_people=10, rows_each=5, seed=1)
    assert X.shape == (50, 4) and y.shape == (50,) and groups.shape == (50,)
    assert len(set(groups.tolist())) == 10
    # Every row belonging to one person carries that person's single label.
    for person in range(10):
        member_labels = set(y[groups == person].tolist())
        assert len(member_labels) == 1


def test_the_regime_series_really_changes_its_rule():
    X, y = s.regime_series(length=600, n_regimes=3, seed=2)
    assert X.shape == (600, 4) and y.shape == (600,)
    # Each regime is internally consistent but the regimes disagree: a
    # linear model fitted on one block should do worse on another.
    from sklearn.linear_model import LogisticRegression

    first = LogisticRegression(max_iter=1000).fit(X[:200], y[:200])
    own = s.accuracy(y[:200], first.predict(X[:200]))
    other = s.accuracy(y[400:], first.predict(X[400:]))
    assert own > 0.9
    assert other < own - 0.2


def test_the_gated_test_set_counts_and_refuses():
    class AlwaysZero:
        def predict(self, X):
            return np.zeros(len(X), dtype=int)

    y = np.array([0, 0, 0, 1])
    gate = s.GatedTestSet(np.zeros((4, 2)), y)
    assert gate.evaluate(AlwaysZero()) == 0.75
    with pytest.raises(s.TestSetTouchedTwice):
        gate.evaluate(AlwaysZero())
    # A fresh gate is a fresh budget; the class holds no global state.
    assert s.GatedTestSet(np.zeros((4, 2)), y).evaluate(AlwaysZero()) == 0.75

#!/usr/bin/env bash
# Day 144 lab harness: "Three Sets, and Why"
#
# 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}"

# Clear caches at the START so the final cleanliness check measures what
# THIS run left behind, not what a previous `pytest starter -q` left.
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
from sklearn.linear_model import LogisticRegression

import splits_lib as s

errors = []


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


# 1. Selection bias
rows = s.selection_bias_curve([1, 2, 5, 10, 25, 50, 100, 500, 1000])
expect(
    "selection bias curve",
    rows,
    [
        (1, 0.4984, 0.5011, -0.0028),
        (2, 0.5115, 0.4992, 0.0123),
        (5, 0.5256, 0.5005, 0.0251),
        (10, 0.5331, 0.4999, 0.0332),
        (25, 0.5436, 0.5009, 0.0427),
        (50, 0.5508, 0.5014, 0.0493),
        (100, 0.5567, 0.4992, 0.0575),
        (500, 0.5682, 0.4978, 0.0704),
        (1000, 0.572, 0.4992, 0.0728),
    ],
)
validation = [v for _k, v, _t, _o in rows]
if not all(a < b for a, b in zip(validation, validation[1:])):
    errors.append("the validation column was not strictly increasing in K")
test_column = [t for _k, _v, t, _o in rows]
if max(test_column) - min(test_column) >= 0.005:
    errors.append("the test column moved, which it must not")
for score in test_column:
    if abs(score - 0.5) >= 0.005:
        errors.append(f"a test score drifted from chance: {score}")

# 1b. The optimism equals the expected maximum
standard_error = s.proportion_standard_error(0.5, 500)
expect("standard error at n=500", round(standard_error, 4), 0.0224)
for k, _v, _t, optimism in rows[1:]:
    in_errors = optimism / standard_error
    simulated = s.expected_max_of_normals(k)
    if abs(in_errors - simulated) >= 0.2:
        errors.append(f"K={k}: optimism {in_errors:.2f} SEs vs expected max {simulated:.2f}")
    if s.sqrt_two_log_k(k) <= simulated:
        errors.append(f"K={k}: sqrt(2 ln K) did not overestimate the expected maximum")
expect("expected max at K=100", round(s.expected_max_of_normals(100), 2), 2.5)
expect("sqrt(2 ln 100)", round(s.sqrt_two_log_k(100), 2), 3.03)

# 2. Stratification
X, y = s.rare_class_dataset()
expect("population positive rate", round(float(y.mean()), 4), 0.05)
random_rates, stratified_rates, empty = s.split_positive_rates(X, y)
expect(
    "random split spread",
    s.spread(random_rates),
    {"mean": 0.0504, "sd": 0.0265, "min": 0.0, "max": 0.16},
)
expect(
    "stratified split spread",
    s.spread(stratified_rates),
    {"mean": 0.05, "sd": 0.01, "min": 0.04, "max": 0.06},
)
expect("random splits with no positives", empty, 21)

# 3. Groups
Xg, yg, groups = s.grouped_dataset()
rowwise, group_aware = s.rowwise_vs_group_split(Xg, yg, groups)
expect("row-wise split", round(rowwise, 4), 0.976)
expect("group-aware split", round(group_aware, 4), 0.4112)
expect("accuracy invented by ignoring groups", round(rowwise - group_aware, 4), 0.5648)
expect("people in both halves", s.groups_shared_between_halves(groups), 50)
if not (group_aware < 0.5 < rowwise):
    errors.append("the group-aware score was not below chance while the row-wise one was above it")

# 4. Time
temporal = s.temporal_inflation_over_constructions()
inflation = [gap for _s, _sh, _ch, _b, gap in temporal]
expect("constructions", len(temporal), 20)
expect("constructions where shuffling won", sum(1 for g in inflation if g > 0), 20)
expect("shuffled mean", round(float(np.mean([r[1] for r in temporal])), 4), 0.5961)
expect("chronological mean", round(float(np.mean([r[2] for r in temporal])), 4), 0.5233)
expect("baseline mean", round(float(np.mean([r[3] for r in temporal])), 4), 0.5235)
expect("inflation mean", round(float(np.mean(inflation)), 4), 0.0728)
expect("inflation sd", round(float(np.std(inflation)), 4), 0.0596)
expect("smallest inflation", round(min(inflation), 4), 0.016)
expect("largest inflation", round(max(inflation), 4), 0.2557)
expect("ratio of largest to smallest", round(max(inflation) / min(inflation), 1), 16.0)

# 5. Holdout versus cross-validation
Xw, yw = s.weak_signal_dataset()
holdout, cross = s.holdout_vs_cross_validation(Xw, yw)
expect(
    "holdout spread",
    s.spread(holdout),
    {"mean": 0.7519, "sd": 0.0381, "min": 0.66, "max": 0.85},
)
expect(
    "cross-validated spread",
    s.spread(cross),
    {"mean": 0.7546, "sd": 0.0061, "min": 0.7375, "max": 0.77},
)
expect("holdout range", round(max(holdout) - min(holdout), 4), 0.19)
expect("steadiness ratio", round(float(np.std(holdout) / np.std(cross)), 4), 6.2344)

# 6. Test-set size
expect(
    "test size table",
    s.test_size_table([50, 100, 200, 500, 1000, 5000]),
    [
        (50, 0.0505, 0.0505, 0.099),
        (100, 0.0357, 0.0357, 0.07),
        (200, 0.0252, 0.0254, 0.0495),
        (500, 0.016, 0.016, 0.0313),
        (1000, 0.0113, 0.0112, 0.0221),
        (5000, 0.005, 0.0051, 0.0099),
    ],
)
expect("rows for +/-0.02", s.rows_needed_for_precision(0.85, 0.02), 1225)
expect("rows for +/-0.01", s.rows_needed_for_precision(0.85, 0.01), 4899)

# 7. The gate
model = LogisticRegression(max_iter=1000).fit(Xw, yw)
gate = s.GatedTestSet(Xw, yw)
expect("first evaluation", round(gate.evaluate(model), 4), 0.7575)
expect("evaluation counter", gate.evaluations, 1)
try:
    gate.evaluate(model)
except s.TestSetTouchedTwice as exc:
    if "validation score" not in str(exc):
        errors.append(f"the gate's message did not explain itself: {exc}")
    if gate.evaluations != 1:
        errors.append("the counter advanced on a refused evaluation")
else:
    errors.append("the gate permitted a second evaluation, which it must not")

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-7 reproduced directly against splits_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}/d144-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_splits_claims.py" <<'PYEOF'
import sys
path = sys.argv[1]
text = open(path).read()
needle = "assert round(rowwise - group_aware, 4) == 0.5648"
replacement = "assert round(rowwise - group_aware, 4) == 0.0"
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_03_splitting_rows_when_the_unit_is_a_person_invents_fifty_six_points"; then
  ok "breaking exercise 3'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 direction of every split result holds beyond the quoted seed"
DIRECTION=$("$PYTHON" - <<'PYEOF'
import sys
sys.path.insert(0, "examples")
import numpy as np
import splits_lib as s

problems = []

# Group leakage is not a property of one dataset seed.
for seed in range(5):
    X, y, groups = s.grouped_dataset(n_people=40, rows_each=15, seed=seed)
    rowwise, grouped = s.rowwise_vs_group_split(X, y, groups, splits=5)
    if rowwise <= grouped:
        problems.append(f"seed {seed}: row-wise {rowwise:.4f} did not beat group-aware {grouped:.4f}")

# Selection optimism is not a property of one replication count.
short = s.selection_bias_curve([1, 100], replications=120)
if short[1][1] <= short[0][1]:
    problems.append("selection optimism vanished at a different replication count")

# Stratification always narrows the spread.
for seed in (144, 145, 146):
    X, y = s.rare_class_dataset(seed=seed)
    r, st, _empty = s.split_positive_rates(X, y, splits=200)
    if s.spread(st)["sd"] >= s.spread(r)["sd"]:
        problems.append(f"seed {seed}: stratifying did not narrow the spread")

if problems:
    for p in problems:
        print("ERROR:", p)
else:
    print("every direction held")
PYEOF
)
if [ "$DIRECTION" = "every direction held" ]; then
  ok "group leakage, selection optimism and stratification hold at seeds the lesson does not quote"
else
  fail "a direction failed beyond the quoted seed"
  echo "$DIRECTION" | 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