Hands-on lab — Day 154: A Complete Regression Project

"""One regression project, run properly, once.

Days 148-153 each isolated one discipline: the one-predictor model and its
four assumptions, the loss as a choice, multicollinearity, ridge and lasso,
the metrics that can be gamed or inverted, and OLS built from scratch. This
module spends every one of those disciplines on a single real dataset and
produces one defensible verdict, with residual diagnostics as the
centrepiece no other day owns.

Frame, baseline, split, pipeline, cross-validate, select, ONE test
evaluation, residual diagnostics, a fairness check, prediction intervals,
an honest interval on the margin. Nothing here is taught for the first
time; everything here is used.
"""

from __future__ import annotations

import numpy as np

from sklearn.datasets import load_diabetes
from sklearn.dummy import DummyRegressor
from sklearn.linear_model import Lasso, LinearRegression, Ridge
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from sklearn.model_selection import KFold, cross_val_predict, cross_val_score, train_test_split
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

# --------------------------------------------------------------------------
# 1. The dataset -- the only bundled regression set that needs no download
# --------------------------------------------------------------------------


def load_dataset():
    """The Wisconsin/diabetes progression set, in raw measurement units.

    ``load_diabetes`` is the only regression dataset scikit-learn bundles
    offline; ``fetch_california_housing`` downloads and is not used here.
    ``scaled=False`` keeps the ten features in their original units (age in
    years, bmi, average blood pressure, six serum measures) so a coefficient
    would still mean something if this project stopped to interpret one --
    Day 148's point, carried forward. The target itself is a composite
    disease-progression score with no physical unit; it is not measured in
    anything, and this project never pretends otherwise.
    """
    d = load_diabetes(scaled=False)
    return d.data, d.target, list(d.feature_names)


# --------------------------------------------------------------------------
# 2. The frame and the baseline -- before any model
# --------------------------------------------------------------------------


def baseline_metrics(x_train, y_train, x_test, y_test):
    """The mean-predictor baseline: RMSE and R^2, computed before any model.

    Day 141's rule, restated for regression: a score is not evidence until
    you know what it beats. Predicting the training mean for every row is
    the simplest possible non-model.
    """
    dummy = DummyRegressor(strategy="mean").fit(x_train, y_train)
    pred = dummy.predict(x_test)
    rmse = float(np.sqrt(mean_squared_error(y_test, pred)))
    r2 = float(r2_score(y_test, pred))
    return round(rmse, 4), round(r2, 4)


# --------------------------------------------------------------------------
# 3. The split -- 442 rows is small, so the test set stays small too
# --------------------------------------------------------------------------


def split_once(X, y, seed: int = 0, test_size: float = 0.25):
    """One split. The test half is touched once, later, for scoring.

    A 25 percent test set on 442 rows is about 110 rows -- small enough
    that every interval in this project is wide, and that is reported
    honestly rather than smoothed over.
    """
    return train_test_split(X, y, test_size=test_size, random_state=seed)


# --------------------------------------------------------------------------
# 4. The candidate pipelines -- ridge, lasso and plain OLS, K counted
# --------------------------------------------------------------------------

_ALPHAS = [0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30, 100]


def candidate_configs():
    """23 candidate pipelines: 11 ridge, 11 lasso, 1 plain OLS.

    Every candidate is a scikit-learn ``Pipeline`` with a ``StandardScaler``
    ahead of the estimator, so cross-validation refits the scaler on each
    fold's training rows only -- Day 143's stage-ordering rule, enforced by
    the estimator's own contract. Returns ``(family, hyperparameter,
    make_pipeline)`` where ``make_pipeline`` is a zero-argument callable
    returning a fresh, unfitted ``Pipeline``.
    """
    configs = []
    for a in _ALPHAS:
        configs.append(
            ("ridge", a, lambda a=a: Pipeline([("scale", StandardScaler()), ("clf", Ridge(alpha=a))]))
        )
    for a in _ALPHAS:
        configs.append(
            (
                "lasso",
                a,
                lambda a=a: Pipeline(
                    [("scale", StandardScaler()), ("clf", Lasso(alpha=a, max_iter=20000))]
                ),
            )
        )
    configs.append(
        ("ols", 0.0, lambda: Pipeline([("scale", StandardScaler()), ("clf", LinearRegression())]))
    )
    return configs


def candidate_count() -> int:
    """K, the number of configurations this project actually tries."""
    return len(candidate_configs())


# --------------------------------------------------------------------------
# 5. Cross-validate, then select -- honest spending of the train rows
# --------------------------------------------------------------------------


def cross_validate_configs(x_train, y_train, seed: int = 0, folds: int = 5):
    """5-fold CV RMSE for every candidate, on train rows only.

    Returns rows of ``(family, hyperparameter, cv_rmse, cv_std)``, sorted
    best (lowest RMSE) first. RMSE is chosen as the selection metric because
    it is stated as the metric BEFORE any model is fitted -- Day 152's whole
    argument that the metric is a choice, made here in the units the target
    is actually in (unitless composite-score points, not mg/dL).
    """
    splitter = KFold(n_splits=folds, shuffle=True, random_state=seed)
    rows = []
    for family, param, make in candidate_configs():
        scores = cross_val_score(make(), x_train, y_train, cv=splitter, scoring="neg_root_mean_squared_error")
        rows.append((family, param, round(float(-scores.mean()), 4), round(float(scores.std()), 4)))
    rows.sort(key=lambda r: r[2])
    return rows


def select_best(x_train, y_train, seed: int = 0, folds: int = 5):
    """Fit the winner of the sweep on the full training set.

    Returns ``(family, hyperparameter, cv_rmse, fitted_pipeline)``.
    """
    rows = cross_validate_configs(x_train, y_train, seed=seed, folds=folds)
    winner_family, winner_param, winner_cv, _sd = rows[0]
    for family, param, make_fn in candidate_configs():
        if family == winner_family and param == winner_param:
            fitted = make_fn().fit(x_train, y_train)
            return winner_family, winner_param, winner_cv, fitted
    raise RuntimeError("winning configuration vanished between sweep and refit")


# --------------------------------------------------------------------------
# 6. The test set -- one look, enforced mechanically
# --------------------------------------------------------------------------


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


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

    Day 144's discipline, reused unchanged: the counter does not advance on
    a refused attempt, so a caller that never succeeds cannot drain the
    budget by retrying. ``evaluate`` returns ``(rmse, r2, mae)`` -- the one
    look this project spends. Reading predictions afterward for residual
    diagnostics is inspection, not a second selection, exactly as Day 147's
    confusion matrix read predictions after its own one evaluation.
    """

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

    def evaluate(self, model):
        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
        pred = model.predict(self._X)
        rmse = float(np.sqrt(mean_squared_error(self._y, pred)))
        r2 = float(r2_score(self._y, pred))
        mae = float(mean_absolute_error(self._y, pred))
        return round(rmse, 4), round(r2, 4), round(mae, 4)


# --------------------------------------------------------------------------
# 7. The verdict -- a bootstrap interval around the margin
# --------------------------------------------------------------------------


def margin_bootstrap_interval(y_test, pred_baseline, pred_model, n_boot: int = 2000, seed: int = 0):
    """A 95 percent interval around the margin: baseline RMSE minus model RMSE.

    Resamples the test rows with replacement ``n_boot`` times and recomputes
    both RMSEs on each resample, so the interval reflects how much the
    margin itself could plausibly move at this test-set size -- the same
    question Day 144's accuracy interval asked, answered here for RMSE
    without an analytic formula, because RMSE's sampling distribution has
    no equally simple closed form.
    """
    rng = np.random.default_rng(seed)
    y_test = np.asarray(y_test)
    pred_baseline = np.asarray(pred_baseline)
    pred_model = np.asarray(pred_model)
    n = len(y_test)
    margins = np.empty(n_boot)
    for i in range(n_boot):
        idx = rng.integers(0, n, size=n)
        rmse_base = np.sqrt(mean_squared_error(y_test[idx], pred_baseline[idx]))
        rmse_model = np.sqrt(mean_squared_error(y_test[idx], pred_model[idx]))
        margins[i] = rmse_base - rmse_model
    lower = round(float(np.percentile(margins, 2.5)), 4)
    upper = round(float(np.percentile(margins, 97.5)), 4)
    return lower, upper


def margin_distinguishable(lower: float, upper: float) -> bool:
    """Whether the bootstrap interval around the margin excludes zero.

    If the interval spans zero, the honest verdict is "cannot distinguish
    this model from the baseline at this test-set size" -- Day 144's
    cautionary case, now checked for a real margin instead of assumed away.
    """
    return lower > 0.0


# --------------------------------------------------------------------------
# 8. Residual diagnostics -- the centrepiece this day owns
# --------------------------------------------------------------------------


def _normal_ppf(p):
    """The inverse standard-normal CDF, by Acklam's rational approximation.

    Built from scratch because scipy is not one of this lab's three pinned
    dependencies -- Day 153's habit, applied again: when a tool is not
    available, build the minimal piece you need rather than reach past the
    pins. Accurate to about 1.15e-9 across the open interval (0, 1); good
    enough for a Q-Q correlation on 111 points. Public-domain algorithm,
    reference in sources.yml.
    """
    a = [-3.969683028665376e01, 2.209460984245205e02, -2.759285104469687e02,
         1.383577518672690e02, -3.066479806614716e01, 2.506628277459239e00]
    b = [-5.447609879822406e01, 1.615858368580409e02, -1.556989798598866e02,
         6.680131188771972e01, -1.328068155288572e01]
    c = [-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e00,
         -2.549732539343734e00, 4.374664141464968e00, 2.938163982698783e00]
    d = [7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e00,
         3.754408661907416e00]
    p_low = 0.02425
    p = np.asarray(p, dtype=float)
    out = np.empty_like(p)
    low = p < p_low
    high = p > (1 - p_low)
    mid = ~(low | high)
    ql = np.sqrt(-2 * np.log(p[low]))
    out[low] = (((((c[0] * ql + c[1]) * ql + c[2]) * ql + c[3]) * ql + c[4]) * ql + c[5]) / (
        (((d[0] * ql + d[1]) * ql + d[2]) * ql + d[3]) * ql + 1
    )
    q = p[mid] - 0.5
    r = q * q
    out[mid] = (((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * q / (
        (((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1)
    )
    qh = np.sqrt(-2 * np.log(1 - p[high]))
    out[high] = -(((((c[0] * qh + c[1]) * qh + c[2]) * qh + c[3]) * qh + c[4]) * qh + c[5]) / (
        (((d[0] * qh + d[1]) * qh + d[2]) * qh + d[3]) * qh + 1
    )
    return out


def residual_summary(y_test, pred_test):
    """Mean and standard deviation of the test residuals, y minus prediction."""
    resid = np.asarray(y_test, dtype=float) - np.asarray(pred_test, dtype=float)
    return round(float(resid.mean()), 4), round(float(resid.std()), 4)


def heteroscedasticity_signal(pred_test, y_test):
    """Correlation between the fitted value and the absolute residual.

    A value near zero means the spread of errors does not grow with the
    predicted value -- the "fanning out" pattern a residual-vs-fitted plot
    is built to catch. This is the from-scratch numeric version of reading
    that plot.
    """
    pred_test = np.asarray(pred_test, dtype=float)
    resid = np.asarray(y_test, dtype=float) - pred_test
    corr = np.corrcoef(pred_test, np.abs(resid))[0, 1]
    return round(float(corr), 4)


def curvature_signal(pred_test, y_test):
    """Correlation between the squared fitted value and the signed residual.

    A value far from zero suggests the model is missing a systematic curve
    -- residuals that trend up then down (or the reverse) as the fitted
    value rises, rather than scattering evenly around zero.
    """
    pred_test = np.asarray(pred_test, dtype=float)
    resid = np.asarray(y_test, dtype=float) - pred_test
    corr = np.corrcoef(pred_test**2, resid)[0, 1]
    return round(float(corr), 4)


def normal_probability_correlation(y_test, pred_test):
    """A from-scratch normal-probability (Q-Q) check, as one correlation.

    Standardises the residuals, sorts them, and correlates them against the
    theoretical normal quantiles a perfectly Gaussian set of residuals
    would produce. 1.0 is a perfectly straight Q-Q line; values noticeably
    below 1.0 flag departures from normality that a plotted Q-Q line would
    show as curvature at the tails.
    """
    resid = np.asarray(y_test, dtype=float) - np.asarray(pred_test, dtype=float)
    n = len(resid)
    std_resid = (resid - resid.mean()) / resid.std()
    sorted_resid = np.sort(std_resid)
    probs = (np.arange(1, n + 1) - 0.5) / n
    theoretical = _normal_ppf(probs)
    corr = np.corrcoef(theoretical, sorted_resid)[0, 1]
    return round(float(corr), 4)


def largest_residuals(y_test, pred_test, n: int = 5):
    """The n largest-magnitude residuals, inspected individually.

    Returns rows of ``(test_row_index, true_value, predicted_value,
    residual)``, sorted by |residual| descending. A confusion matrix reads
    the specific mistakes a classifier makes; this is that discipline's
    regression counterpart.
    """
    y_test = np.asarray(y_test, dtype=float)
    pred_test = np.asarray(pred_test, dtype=float)
    resid = y_test - pred_test
    order = np.argsort(-np.abs(resid))[:n]
    return [
        (int(i), round(float(y_test[i]), 4), round(float(pred_test[i]), 4), round(float(resid[i]), 4))
        for i in order
    ]


# --------------------------------------------------------------------------
# 9. Is the model worse for high-value targets? Measure it.
# --------------------------------------------------------------------------


def error_by_target_level(y_test, pred_test):
    """RMSE on the below-median half of test targets against the above-median half.

    Returns ``(rmse_low, rmse_high, ratio)`` where ``ratio`` is
    ``rmse_high / rmse_low``. On a disease-progression score this is a
    fairness-relevant question, not only a statistical one: are errors
    worse for patients whose true progression is more severe?
    """
    y_test = np.asarray(y_test, dtype=float)
    pred_test = np.asarray(pred_test, dtype=float)
    median_y = float(np.median(y_test))
    low_mask = y_test <= median_y
    high_mask = ~low_mask
    rmse_low = float(np.sqrt(mean_squared_error(y_test[low_mask], pred_test[low_mask])))
    rmse_high = float(np.sqrt(mean_squared_error(y_test[high_mask], pred_test[high_mask])))
    return round(rmse_low, 4), round(rmse_high, 4), round(rmse_high / rmse_low, 4)


# --------------------------------------------------------------------------
# 10. The leaky version -- selecting by peeking at the test set
# --------------------------------------------------------------------------


def leaky_selection_test_rmse(x_train, y_train, x_test, y_test):
    """Select the winner by fitting every candidate and scoring it on TEST.

    Returns the winning (lowest) test RMSE -- K looks disguised as one,
    Day 147's mistake, reconstructed for regression. Lower RMSE is better,
    so the leaky search can only match or beat the honestly selected
    model's own test RMSE, never lose to it.
    """
    best_rmse = None
    for _family, _param, make in candidate_configs():
        pipe = make().fit(x_train, y_train)
        pred = pipe.predict(x_test)
        rmse = float(np.sqrt(mean_squared_error(y_test, pred)))
        if best_rmse is None or rmse < best_rmse:
            best_rmse = rmse
    return round(best_rmse, 4)


def leaky_vs_honest_over_seeds(X, y, seeds=range(20), folds: int = 5):
    """The gap between peeking at the test set and looking at it once.

    Returns rows of ``(seed, honest_rmse, leaky_rmse, gap)`` where
    ``gap = honest_rmse - leaky_rmse``. Because lower RMSE is better, a
    positive gap means the leak reported a lower (better-looking) error
    than the honest evaluation -- the leak can only help the reported
    number, never hurt it.
    """
    rows = []
    for seed in seeds:
        x_train, x_test, y_train, y_test = split_once(X, y, seed=seed)
        _family, _param, _cv, fitted = select_best(x_train, y_train, seed=seed, folds=folds)
        honest_rmse = float(np.sqrt(mean_squared_error(y_test, fitted.predict(x_test))))
        leaky_rmse = leaky_selection_test_rmse(x_train, y_train, x_test, y_test)
        rows.append((seed, round(honest_rmse, 4), leaky_rmse, round(honest_rmse - leaky_rmse, 4)))
    return rows


# --------------------------------------------------------------------------
# 11. Prediction intervals -- not just a point, and measured coverage
# --------------------------------------------------------------------------


def prediction_interval_coverage(x_train, y_train, x_test, y_test, fitted, seed: int = 0, folds: int = 5):
    """A constant-width 95 percent prediction interval, and its realised coverage.

    The half-width comes from the standard deviation of out-of-fold
    residuals on the TRAINING rows only (``cross_val_predict``), never from
    the test residuals themselves -- using the test residuals to size the
    test interval would be circular. Returns ``(half_width, coverage)``
    where ``coverage`` is the fraction of test targets that actually fall
    inside ``prediction +/- half_width``, measured against the nominal 0.95.
    """
    splitter = KFold(n_splits=folds, shuffle=True, random_state=seed)
    oof_pred = cross_val_predict(fitted, x_train, y_train, cv=splitter)
    oof_resid = np.asarray(y_train, dtype=float) - oof_pred
    half_width = round(float(1.96 * oof_resid.std()), 4)
    pred_test = fitted.predict(x_test)
    y_test = np.asarray(y_test, dtype=float)
    lower = pred_test - half_width
    upper = pred_test + half_width
    coverage = round(float(np.mean((y_test >= lower) & (y_test <= upper))), 4)
    return half_width, coverage

#!/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.dummy import DummyRegressor  # noqa: E402

import regression_lib as r  # noqa: E402


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


def main() -> None:
    print("Day 154 -- a complete regression project, measured")
    print("=" * 51)

    rule("1. The dataset: the only bundled regression set")
    X, y, names = r.load_dataset()
    print(f"  shape: {X.shape}  features: {names}")
    print(f"  y range: [{float(y.min()):.1f}, {float(y.max()):.1f}]  mean: {float(y.mean()):.4f}")
    print("  the target is a composite disease-progression score with no physical unit")

    x_train, x_test, y_train, y_test = r.split_once(X, y, seed=0)
    print(f"  split: train={x_train.shape[0]}  test={x_test.shape[0]}  (seed 0, 25 percent test)")

    rule("2. The baseline, before any model")
    base_rmse, base_r2 = r.baseline_metrics(x_train, y_train, x_test, y_test)
    print(f"  mean-predictor baseline: RMSE {base_rmse:.4f}  R2 {base_r2:.4f}")

    rule("3. The sweep: cross-validate every candidate on train rows only")
    k = r.candidate_count()
    print(f"  K = {k} candidate pipelines: 11 ridge, 11 lasso, 1 plain OLS")
    family, param, cv_rmse, fitted = r.select_best(x_train, y_train, seed=0)
    print(f"  winner: {family} (alpha={param})  5-fold CV RMSE = {cv_rmse:.4f}")

    rule("4. ONE test evaluation")
    gate = r.GatedTestSet(x_test, y_test)
    test_rmse, test_r2, test_mae = gate.evaluate(fitted)
    print(f"  test RMSE: {test_rmse:.4f}  R2: {test_r2:.4f}  MAE: {test_mae:.4f}")
    try:
        gate.evaluate(fitted)
        print("  second evaluation: NO ERROR RAISED")
    except r.TestSetTouchedTwice as exc:
        print(f"  second evaluation : {type(exc).__name__}")
        print(f"    {exc}")

    rule("5. The margin, with a bootstrap interval")
    baseline_model = DummyRegressor(strategy="mean").fit(x_train, y_train)
    pred_baseline = baseline_model.predict(x_test)
    pred_model = fitted.predict(x_test)
    margin = round(base_rmse - test_rmse, 4)
    lower, upper = r.margin_bootstrap_interval(y_test, pred_baseline, pred_model, seed=0)
    print(f"  margin (baseline RMSE - model RMSE): {margin:+.4f}")
    print(f"  95 percent bootstrap interval on the margin: [{lower:.4f}, {upper:.4f}]")
    print(f"  distinguishable from baseline at this test-set size: {r.margin_distinguishable(lower, upper)}")

    rule("6. Residual diagnostics -- the centrepiece")
    resid_mean, resid_std = r.residual_summary(y_test, pred_model)
    het = r.heteroscedasticity_signal(pred_model, y_test)
    curv = r.curvature_signal(pred_model, y_test)
    qq = r.normal_probability_correlation(y_test, pred_model)
    print(f"  residuals: mean {resid_mean:+.4f}  sd {resid_std:.4f}")
    print(f"  heteroscedasticity signal, corr(fitted, |residual|): {het:+.4f}")
    print(f"  curvature signal, corr(fitted^2, residual): {curv:+.4f}")
    print(f"  normal-probability (Q-Q) correlation: {qq:.4f}")
    print("  largest residuals (row, true, predicted, residual):")
    for row_idx, true_val, pred_val, resid in r.largest_residuals(y_test, pred_model, n=5):
        print(f"    row {row_idx:3d}: true={true_val:7.1f}  pred={pred_val:7.1f}  residual={resid:+8.1f}")

    rule("7. Is the model worse for high-value targets?")
    rmse_low, rmse_high, ratio = r.error_by_target_level(y_test, pred_model)
    print(f"  RMSE on below-median targets: {rmse_low:.4f}")
    print(f"  RMSE on above-median targets: {rmse_high:.4f}")
    print(f"  ratio (high / low): {ratio:.4f}")

    rule("8. The leaky version: selecting by peeking at the test set")
    leaky_rmse = r.leaky_selection_test_rmse(x_train, y_train, x_test, y_test)
    print(f"  honest (select on CV, look once): {test_rmse:.4f}")
    print(f"  leaky (best of {k} scored directly on test): {leaky_rmse:.4f}")
    print(f"  gap (honest - leaky, positive means the leak looked better): {round(test_rmse - leaky_rmse, 4):+.4f}")

    rule("8b. The leaky gap, over 20 seeds")
    rows = r.leaky_vs_honest_over_seeds(X, y, seeds=range(20))
    gaps = np.array([row[3] for row in rows])
    print("   seed   honest    leaky      gap")
    for seed, honest_s, leaky_s, gap_s in rows:
        print(f"  {seed:5d}   {honest_s:.4f}   {leaky_s:.4f}   {gap_s:+.4f}")
    print(f"  mean gap: {float(gaps.mean()):+.4f}  sd {float(gaps.std()):.4f}  min {float(gaps.min()):+.4f}  max {float(gaps.max()):+.4f}")
    print(f"  fraction of seeds where the leak was non-negative: {float((gaps >= 0).mean()):.4f}")

    rule("9. Prediction intervals, and their realised coverage")
    half_width, coverage = r.prediction_interval_coverage(x_train, y_train, x_test, y_test, fitted, seed=0)
    print(f"  95 percent prediction interval half-width (from TRAIN out-of-fold residuals): +/-{half_width:.4f}")
    print(f"  realised coverage on the 111 test rows: {coverage:.4f}  (nominal: 0.9500)")

    rule("10. What the whole thing costs")
    print("  wall-clock cost is machine-dependent and not reproduced here byte for byte;")
    print("  see metadata.yml and expected-output/FIELDS.md for the captured timing")


if __name__ == "__main__":
    main()

"""The reference solutions: one regression project, run properly, once.

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 inspect

import numpy as np
import pytest
from sklearn.datasets import fetch_california_housing, load_diabetes
from sklearn.dummy import DummyRegressor

import regression_lib as r


@pytest.fixture(scope="module")
def dataset():
    return r.load_dataset()


@pytest.fixture(scope="module")
def split(dataset):
    X, y, _names = dataset
    return r.split_once(X, y, seed=0)


# --- 1. The dataset, and why it is the only one used ------------------------


def test_01_the_dataset_is_the_only_bundled_regression_set(dataset):
    X, y, names = dataset
    assert X.shape == (442, 10)
    assert names == ["age", "sex", "bmi", "bp", "s1", "s2", "s3", "s4", "s5", "s6"]
    assert round(float(y.min()), 4) == 25.0
    assert round(float(y.max()), 4) == 346.0
    assert round(float(y.mean()), 4) == 152.1335
    # fetch_california_housing downloads by default -- forbidden by this
    # lab's offline rule -- while load_diabetes never touches the network.
    sig = inspect.signature(fetch_california_housing)
    assert sig.parameters["download_if_missing"].default is True


def test_01b_raw_units_are_not_the_scikit_learn_default(dataset):
    X, y, names = dataset
    # This lab's loader asks for scaled=False. Compare against scikit-learn's
    # own default (scaled=True) to see why: the default is mean-centred and
    # variance-scaled to tiny floats, not the years and mg/dL a reader could
    # sanity-check a coefficient against.
    default = load_diabetes(scaled=True)
    age_index = names.index("age")
    assert X[:, age_index].min() == 19.0
    assert X[:, age_index].max() == 79.0
    assert -0.2 < default.data[:, age_index].min() < 0
    assert 0 < default.data[:, age_index].max() < 0.2


# --- 2. The frame and the baseline -----------------------------------------


def test_02_the_baseline_before_any_model(split):
    x_train, x_test, y_train, y_test = split
    rmse, r2 = r.baseline_metrics(x_train, y_train, x_test, y_test)
    assert rmse == 70.4637
    assert r2 == -0.0001
    # Every model below has to beat 70.4637 RMSE to be worth building.


# --- 3. The split -------------------------------------------------------


def test_03_the_split_holds_the_test_rows_back(dataset):
    X, y, _names = dataset
    x_train, x_test, y_train, y_test = r.split_once(X, y, seed=0)
    assert x_train.shape == (331, 10)
    assert x_test.shape == (111, 10)
    # 111 test rows on a 442-row dataset: small, and every interval below
    # is wide because of it -- that is the point, not a flaw to fix.


# --- 4. The sweep --------------------------------------------------------


def test_04_the_sweep_counts_twenty_three_candidate_pipelines():
    assert r.candidate_count() == 23
    families = [family for family, _param, _make in r.candidate_configs()]
    assert families.count("ridge") == 11
    assert families.count("lasso") == 11
    assert families.count("ols") == 1


# --- 5. Cross-validate, then select --------------------------------------


def test_05_cross_validation_selects_the_winner_on_train_rows_only(split):
    x_train, _x_test, y_train, _y_test = split
    family, param, cv_rmse, fitted = r.select_best(x_train, y_train, seed=0)
    assert (family, param) == ("lasso", 1)
    assert cv_rmse == 53.8958
    assert hasattr(fitted, "predict")
    # The winner was chosen on cross-validated train rows; test has not
    # been touched.


# --- 6. The gate: one look at the test set ----------------------------------


def test_06_the_gate_permits_exactly_one_test_evaluation(split):
    x_train, x_test, y_train, y_test = split
    _family, _param, _cv, fitted = r.select_best(x_train, y_train, seed=0)
    gate = r.GatedTestSet(x_test, y_test)
    assert gate.evaluations == 0
    rmse, r2, mae = gate.evaluate(fitted)
    assert rmse == 56.5566
    assert r2 == 0.3557
    assert mae == 45.2846
    assert gate.evaluations == 1
    with pytest.raises(r.TestSetTouchedTwice) as excinfo:
        gate.evaluate(fitted)
    assert "validation score" in str(excinfo.value)
    assert gate.evaluations == 1


# --- 7. The margin, with a bootstrap interval ------------------------------


def test_07_the_margin_has_a_bootstrap_interval(split):
    x_train, x_test, y_train, y_test = split
    baseline = DummyRegressor(strategy="mean").fit(x_train, y_train)
    _family, _param, _cv, fitted = r.select_best(x_train, y_train, seed=0)
    pred_baseline = baseline.predict(x_test)
    pred_model = fitted.predict(x_test)
    lower, upper = r.margin_bootstrap_interval(y_test, pred_baseline, pred_model, seed=0)
    assert lower == 5.5852
    assert upper == 22.3324
    rmse_base, _r2_base = r.baseline_metrics(x_train, y_train, x_test, y_test)
    rmse_model = float(np.sqrt(np.mean((y_test - pred_model) ** 2)))
    margin = round(rmse_base - rmse_model, 4)
    assert margin == 13.9071
    # The margin (13.9071) clears the interval's lower bound (5.5852): this
    # model IS distinguishable from the baseline at this test-set size.
    assert r.margin_distinguishable(lower, upper) is True


# --- 8. Residual diagnostics -- the centrepiece ----------------------------


def test_08_the_residual_vs_fitted_diagnostic(split):
    x_train, x_test, y_train, y_test = split
    _family, _param, _cv, fitted = r.select_best(x_train, y_train, seed=0)
    pred = fitted.predict(x_test)
    resid_mean, resid_std = r.residual_summary(y_test, pred)
    assert resid_mean == -3.6262
    assert resid_std == 56.4402
    het = r.heteroscedasticity_signal(pred, y_test)
    curv = r.curvature_signal(pred, y_test)
    assert het == 0.2386
    assert curv == -0.1278
    # A modest positive heteroscedasticity signal (errors fan out a little
    # as predictions rise) and a weak curvature signal (no strong missed
    # trend) -- neither dramatic, both worth reporting rather than ignoring.


def test_08b_the_normal_probability_check_and_the_largest_residuals(split):
    x_train, x_test, y_train, y_test = split
    _family, _param, _cv, fitted = r.select_best(x_train, y_train, seed=0)
    pred = fitted.predict(x_test)
    qq = r.normal_probability_correlation(y_test, pred)
    assert qq == 0.9901
    rows = r.largest_residuals(y_test, pred, n=5)
    assert len(rows) == 5
    assert rows[0] == (60, 52.0, 209.3314, -157.3314)
    assert rows[1] == (65, 302.0, 153.8865, 148.1135)
    # 0.9901 is close to a perfectly straight Q-Q line: these residuals are
    # not wildly non-normal, even on real data with only 111 test rows.


# --- 9. Is the model worse for high-value targets? Measure it. -------------


def test_09_error_by_target_level(split):
    x_train, x_test, y_train, y_test = split
    _family, _param, _cv, fitted = r.select_best(x_train, y_train, seed=0)
    pred = fitted.predict(x_test)
    rmse_low, rmse_high, ratio = r.error_by_target_level(y_test, pred)
    assert rmse_low == 55.2464
    assert rmse_high == 57.8601
    assert ratio == 1.0473
    # Only 4.73 percent worse on the more-severe half: not the dramatic
    # fairness problem this exercise sets out to check for, at this seed --
    # and that is exactly why the check has to run rather than be assumed.


# --- 10. The leaky version --------------------------------------------------


def test_10_the_leaky_version_selects_by_peeking_at_the_test_set(split):
    x_train, x_test, y_train, y_test = split
    _family, _param, _cv, fitted = r.select_best(x_train, y_train, seed=0)
    honest_rmse = round(float(np.sqrt(np.mean((y_test - fitted.predict(x_test)) ** 2))), 4)
    leaky_rmse = r.leaky_selection_test_rmse(x_train, y_train, x_test, y_test)
    assert honest_rmse == 56.5566
    assert leaky_rmse == 55.5212
    # Lower RMSE is better: the leak can only match or beat the honest
    # score, never lose to it.
    assert leaky_rmse <= honest_rmse


def test_10b_the_leaky_gap_over_twenty_seeds(dataset):
    X, y, _names = dataset
    rows = r.leaky_vs_honest_over_seeds(X, y, seeds=range(20))
    assert len(rows) == 20
    gaps = np.array([row[3] for row in rows])
    assert round(float(gaps.mean()), 4) == 0.5279
    assert round(float(gaps.std()), 4) == 0.3686
    assert round(float(gaps.min()), 4) == 0.011
    assert round(float(gaps.max()), 4) == 1.1451
    # The leak never once hurt the reported number, across 20 independent
    # seeds: the mechanism, not luck.
    assert (gaps >= 0).all()


# --- 11. Prediction intervals, and their realised coverage -----------------


def test_11_prediction_interval_coverage(split):
    x_train, x_test, y_train, y_test = split
    _family, _param, _cv, fitted = r.select_best(x_train, y_train, seed=0)
    half_width, coverage = r.prediction_interval_coverage(
        x_train, y_train, x_test, y_test, fitted, seed=0
    )
    assert half_width == 105.8797
    assert coverage == 0.9459
    # Nominal is 0.95; measured is 0.9459 -- close, on only 111 test rows.
    assert 0.85 < coverage <= 1.0

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

These five 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 regression_lib as r


def test_the_dataset_loads_offline_and_matches_its_shape():
    X, y, names = r.load_dataset()
    assert X.shape == (442, 10)
    assert y.shape == (442,)
    assert names == ["age", "sex", "bmi", "bp", "s1", "s2", "s3", "s4", "s5", "s6"]
    assert float(y.min()) == 25.0
    assert float(y.max()) == 346.0


def test_candidate_configs_really_are_twenty_three_distinct_pipelines():
    configs = r.candidate_configs()
    assert len(configs) == 23
    seen = set()
    for family, param, make in configs:
        pipe = make()
        assert hasattr(pipe, "fit") and hasattr(pipe, "predict")
        seen.add((family, param))
    # No two configs share a (family, hyperparameter) pair.
    assert len(seen) == 23


def test_the_gated_test_set_counts_and_refuses():
    class AlwaysMean:
        def predict(self, X):
            return np.full(len(X), 150.0)

    y = np.array([100.0, 200.0, 150.0, 150.0])
    gate = r.GatedTestSet(np.zeros((4, 2)), y)
    rmse, r2, mae = gate.evaluate(AlwaysMean())
    assert rmse > 0
    with pytest.raises(r.TestSetTouchedTwice):
        gate.evaluate(AlwaysMean())
    # A fresh gate is a fresh budget; the class holds no global state.
    fresh_rmse, _r2, _mae = r.GatedTestSet(np.zeros((4, 2)), y).evaluate(AlwaysMean())
    assert fresh_rmse == rmse


def test_the_normal_ppf_matches_known_reference_points():
    # Well-known standard-normal quantiles, to a few decimals.
    assert abs(r._normal_ppf(np.array([0.5]))[0] - 0.0) < 1e-6
    assert abs(r._normal_ppf(np.array([0.975]))[0] - 1.959964) < 1e-4
    assert abs(r._normal_ppf(np.array([0.025]))[0] - (-1.959964)) < 1e-4


def test_largest_residuals_returns_them_sorted_by_magnitude():
    y_test = np.array([10.0, 20.0, 30.0, 40.0])
    pred_test = np.array([10.0, 25.0, 10.0, 41.0])
    rows = r.largest_residuals(y_test, pred_test, n=2)
    assert len(rows) == 2
    # Row 2 (residual +20) is the largest in magnitude; row 3 (residual -1) is smallest.
    assert rows[0][0] == 2
    assert rows[0][3] == 20.0

lesson_id: D154
day: 154
kind: guided-build
languages:
  - python
  - bash
setup_commands:
  - cd labs/sections/machine-learning/day-154-a-complete-regression-project
  - 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: 70
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 -> 15 checks, 0 failure(s), exit 0. pytest examples -q -> 19
  passed. pytest starter -q -> 5 passed, 14 skipped (the five machinery checks in
  test_regression_lib.py are solved in both directories; the fourteen exercise stubs
  in starter/test_regression_claims.py are untouched). Everything ran through a real
  lab-local .venv created by the documented setup commands; scikit-learn pulled in
  scipy, joblib and threadpoolctl as its own dependencies, none of which this lab
  imports directly (the Q-Q normal-probability check is built from scratch instead of
  importing scipy). The lab is fully offline after the pip install -- the dataset is
  bundled inside scikit-learn itself (load_diabetes), nothing is downloaded, and
  harness check 10 confirms no URL appears anywhere in starter/ or examples/ source.
  Section 8 of the harness copies examples/ into a mktemp -d scratch directory,
  confirms 19 passed, rewrites the largest-residuals assertion in exercise 8b to a
  value that cannot hold, confirms a non-zero exit naming the failing test, and
  removes the scratch directory. Separately, by hand, the winning-configuration
  assertion in exercise 5 was changed from ('lasso', 1) to ('ridge', 1) in
  examples/test_regression_claims.py and the whole harness re-run: it reported 15
  checks, 2 failure(s) and exited 1 (both the pytest run in section 3 and the
  pytest-free direct reproduction in section 2 caught it); the file was restored and
  the harness returned to 15 checks, 0 failure(s), exit 0. Measured wall-clock time
  for one seed's frame-to-verdict pipeline: 0.1386s. The 20-seed leaky-gap comparison
  in exercise 10b: 2.9559s. The full report_measurements.py run (which performs that
  20-seed sweep and a 2000-draw bootstrap): 3.99s user CPU / 4.13s wall-clock.
  THE DATASET: sklearn.datasets.load_diabetes(scaled=False), the only regression
  dataset bundled inside scikit-learn -- fetch_california_housing downloads by
  default (download_if_missing=True) and is forbidden by this lab's offline rule.
  442 rows, 10 real-valued measurements in raw units (age in years, not a
  mean-centred fraction), target range [25.0, 346.0], mean 152.1335. THE TARGET HAS
  NO PHYSICAL UNIT -- it is a composite disease-progression score, and this lab says
  so plainly rather than inventing one, following Day 152's precedent. MEASURED
  PAIRS, all captured verbatim in expected-output/measured-values.txt. (1) Split
  (seed 0, 25 percent test): 331 train rows, 111 test rows. (2) Mean-predictor
  baseline, computed BEFORE any model: RMSE 70.4637, R2 -0.0001. (3) K = 23 candidate
  pipelines (11 ridge alphas, 11 lasso alphas, 1 plain OLS), 5-fold cross-validation
  on train rows only, scored by RMSE (the metric chosen before modelling, in the
  target's own unitless composite-score points); the winner at seed 0 is
  Lasso(alpha=1) at cv_rmse 53.8958 -- plain OLS scores 54.0926 at the same seed, only
  0.1968 RMSE points behind, so regularisation's benefit on PREDICTION accuracy here
  is real but modest, even though Day 150 already showed multicollinearity makes the
  raw OLS COEFFICIENTS unstable. (4) The ONE test evaluation, enforced by a
  GatedTestSet that raises TestSetTouchedTwice on any further attempt without
  advancing its counter: test RMSE 56.5566, R2 0.3557, MAE 45.2846. (5) The margin
  over baseline (70.4637 - 56.5566 = 13.9071 RMSE points) with a 2000-draw bootstrap
  95 percent interval of [5.5852, 22.3324] (seeded, deterministic under the pinned
  NumPy) -- the interval excludes zero, so the honest verdict is DISTINGUISHABLE, not
  the "cannot distinguish" case a smaller improvement or a smaller test set could have
  forced. (6) RESIDUAL DIAGNOSTICS, the centrepiece: residual mean -3.6262, sd
  56.4402; heteroscedasticity signal (corr of fitted value with absolute residual)
  +0.2386 -- a modest fanning-out, not dramatic; curvature signal (corr of squared
  fitted value with signed residual) -0.1278 -- weak, no strong missed trend; a
  from-scratch Q-Q normal-probability correlation of 0.9901 -- close to a straight
  line, these residuals are close to normal even on real data with only 111 test
  rows. The five largest residuals were inspected individually; the worst is test row
  60, true value 52.0, predicted 209.3314, a residual of -157.3314 points. (7) ERROR
  BY TARGET LEVEL, a fairness-relevant measurement on a disease-progression score:
  RMSE 55.2464 on the below-median half of test targets, 57.8601 on the above-median
  half, a ratio of 1.0473 -- only 4.73 percent worse on the more severe half at this
  seed, not the dramatic asymmetry the check sets out to look for, and that is a real
  finding rather than a disappointing one. (8) The leaky version -- selecting by
  fitting all 23 candidates and scoring each directly on the test set instead of
  selecting on cross-validated train rows: at seed 0 the leaky RMSE is 55.5212 against
  the honest 56.5566, a gap of 1.0354 points in the leak's favour; over 20 seeds the
  mean gap is 0.5279 (sd 0.3686, min 0.011, max 1.1451) and the gap is NEVER negative
  at any seed -- the leak can only make the reported error look as good or better,
  never worse. (9) PREDICTION INTERVALS: a constant 95 percent half-width of 105.8797,
  built from the standard deviation of TRAINING out-of-fold residuals only (never from
  the test residuals, which would be circular); realised coverage on the 111 test
  rows was 0.9459 against the 0.95 nominal rate -- close, and a separate 10-seed check
  (not asserted by the harness, reported here for context) gives a mean coverage of
  0.9558, range 0.9369 to 0.991. FOUR HONESTY CALLS. FIRST: the margin over baseline
  (13.9071 RMSE points) comfortably clears its own bootstrap interval's lower bound
  (5.5852), so this project's verdict is DISTINGUISHABLE -- this is reported as the
  measured outcome, not adjusted to manufacture the "cannot distinguish" case the
  brief flagged as the more valuable one; on this dataset, at this seed, the model
  genuinely beats the baseline by more than sampling noise could explain, and the
  arithmetic that could have said otherwise is run rather than skipped. SECOND: the
  heteroscedasticity and curvature signals (+0.2386 and -0.1278) are real but modest
  -- neither is dramatic, and the lab reports them as such rather than either
  overstating a problem or hiding a real (if mild) fanning-out pattern. THIRD: the
  fairness check by target level came back nearly flat (ratio 1.0473) -- a genuinely
  boring result at this seed, reported plainly rather than searched for a seed where
  it would look more dramatic. FOURTH: the prediction-interval half-width is built
  exclusively from TRAINING out-of-fold residuals, never from the test set's own
  residuals; using test residuals to size an interval evaluated on the same test set
  would be circular, and this lab does not do that even though it would have made
  coverage numbers easier to hit exactly. Harness check 9 re-runs the leaky-gap
  direction and the selection mechanics at seeds this lab does not quote (20-24 and
  seed 41), so the headline directions are confirmed beyond the seeds reported in
  prose.

# 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`. The Q-Q normal-probability check is built
from scratch inside `regression_lib.py` specifically so this lab does not
need scipy.

## Why the versions are pinned exactly

`KFold(shuffle=True, random_state=...)`, `train_test_split`, and the
bootstrap resampling in `margin_bootstrap_interval` all depend on NumPy's
`Generator` bit stream, which NumPy's own documentation states carries no
cross-version compatibility guarantee. Different pinned versions can
legitimately shuffle the same seed into a different order, and every
downstream number -- the winning configuration, its cross-validated RMSE,
the test score, every residual diagnostic -- would move with it.

What does not depend on the pins: the RMSE and R2 formulas, which are
arithmetic; the direction of every structural result -- the leaky
selection score is never worse than the honest one, cross-validation
selects on train rows only, a gated test set refuses a second look; and
structural facts, such as the dataset used here having 442 rows and 10
features. `expected-output/FIELDS.md` separates the two categories in
full.

## 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: the
dataset is bundled inside scikit-learn itself, and nothing else 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 154 lab brief — A Complete Regression Project

Days 148 through 153 each isolated one discipline in a lab built to show
it in isolation: the one-predictor model and its four assumptions, the
loss function as a choice, multicollinearity, ridge and lasso, the
metrics that can be gamed or inverted, and OLS built from scratch.

This lab is not a new discipline. It is all six of them, spent on one
real dataset, in the order a working project actually uses them: frame,
baseline, split, pipeline, cross-validate, select, **one** test
evaluation, residual diagnostics, a fairness check, prediction intervals,
an honest verdict with an interval on the margin.

## The dataset, and why there is only one choice

`sklearn.datasets.load_diabetes` is the only regression dataset that
ships inside scikit-learn and needs no download.
`fetch_california_housing` downloads on first use and is forbidden by
this lab's offline rule. Exercise 1 checks both facts directly.

`load_diabetes(scaled=False)` gives 442 rows of 10 real-valued
measurements — age, sex, bmi, average blood pressure, and six serum
measures — in their **raw units**: age in years, not a mean-centred
fraction. scikit-learn's own default (`scaled=True`) mean-centres and
variance-scales every column to tiny floats; exercise 1b measures the
difference directly.

**The target has no physical unit.** It is a composite disease-
progression score assembled from clinical measurements a year after
baseline — not mg/dL, not a count, not anything with a name. Every RMSE
and MAE in this lab is in those same unitless composite-score points, and
this brief says so plainly rather than inventing a unit that does not
exist.

442 rows is small. A 25 percent test set is 111 rows — small enough that
every interval this lab computes is wide, and that is the point, not a
flaw: an honest verdict on data this size can legitimately come out
"cannot distinguish," and this lab's own margin check (exercise 7)
reports what actually happened rather than assuming an answer.

## What the honest run measures

Twenty-three candidate pipelines — 11 ridge regularisation strengths, 11
lasso regularisation strengths, and 1 plain OLS — are cross-validated
five ways on the training rows only, scored by RMSE. The winner is
`Lasso(alpha=1)`, at a cross-validated RMSE of **53.8958**. Fitted on the
full training set and evaluated **exactly once** against the test rows,
it scores **56.5566** RMSE (R2 0.3557, MAE 45.2846).

## The margin, and whether it survives an interval

The mean-predictor baseline scores 70.4637 RMSE. The winning model's
margin over that baseline is **13.9071** RMSE points. A 2000-draw
bootstrap resample of the 111 test rows gives a 95 percent interval on
that margin of **[5.5852, 22.3324]** — the interval excludes zero, so
this project's honest verdict is that the model IS distinguishable from
the baseline at this test-set size. It could easily have come out the
other way on a smaller improvement or a smaller test set, which is
exactly why the bootstrap runs rather than being assumed.

## Residual diagnostics — the centrepiece of this lab

No other day in this course owns reading the residuals themselves.
Exercise 8 measures whether the errors fan out as predictions rise
(a **heteroscedasticity signal** of 0.2386 — mild, not dramatic), whether
there is a missed curve in the fit (a **curvature signal** of -0.1278 —
weak), and how close the residuals are to normally distributed (a
**Q-Q correlation** of 0.9901, built from scratch with no scipy
dependency, since scipy is not one of this lab's three pinned packages).
Exercise 8b names the five largest individual mistakes: the biggest is a
patient whose true score was 52 and whose prediction was 209.3 — a
residual of -157.3 points.

## Is the model worse for high-value targets? Measured, not assumed.

Exercise 9 splits the test set at the median true value and compares
RMSE on each half: 55.2464 on the below-median half, 57.8601 on the
above-median half — a ratio of **1.0473**. Only 4.73 percent worse on the
more severe half at this seed: not the dramatic fairness problem the
exercise sets out to check for, and that is a real finding, not a
disappointing one — it had to be measured to be known.

## The leak this lab lets you cause on purpose

Exercise 10 rebuilds the mistake Day 147 spent a whole lesson on, for
regression: selecting a model by fitting every candidate and scoring it
**on the test set directly**, keeping whichever scores the lowest RMSE,
instead of selecting on cross-validated train rows and looking at test
once.

At the reported seed the leak reports 55.5212 RMSE against the honest
56.5566 — a gap of 1.0354 points in the leak's favour. Over 20 seeds the
mean gap is **0.5279** (sd 0.3686, min 0.011, max 1.1451), and it is
**never negative** — the leak can only make the reported error look as
good or better than the honest one, never worse. That asymmetry is the
mechanism, not luck.

## Prediction intervals, and their realised coverage

Exercise 11 builds a constant-width 95 percent prediction interval from
the standard deviation of **training** out-of-fold residuals — never from
the test residuals themselves, which would be circular — and checks what
fraction of the 111 test targets actually fall inside it: **0.9459**,
against a 0.95 nominal rate. Close, on 111 rows.

## How to work

1. Build the environment (see the lab `README.md`).
2. Run `.venv/bin/pytest starter -q`. You will see five passes (the
   machinery checks in `test_regression_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. 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 `regression_lib.py`, `test_regression_lib.py` and
`test_regression_claims.py`; pytest aborts on the module-name collision.
Run them separately, always.

## And the rule, made mechanical, again

Exercise 6 wraps the test set in the same `GatedTestSet` pattern Day 144
built and Day 147 reused: exactly one evaluation, `TestSetTouchedTwice`
on the second, and a counter that does not advance on a refused attempt.
Thirteen days in, this is not a new idea — it is the same discipline,
proven on a regression problem instead of a classification one.

"""One regression project, run properly, once.

Days 148-153 each isolated one discipline: the one-predictor model and its
four assumptions, the loss as a choice, multicollinearity, ridge and lasso,
the metrics that can be gamed or inverted, and OLS built from scratch. This
module spends every one of those disciplines on a single real dataset and
produces one defensible verdict, with residual diagnostics as the
centrepiece no other day owns.

Frame, baseline, split, pipeline, cross-validate, select, ONE test
evaluation, residual diagnostics, a fairness check, prediction intervals,
an honest interval on the margin. Nothing here is taught for the first
time; everything here is used.
"""

from __future__ import annotations

import numpy as np

from sklearn.datasets import load_diabetes
from sklearn.dummy import DummyRegressor
from sklearn.linear_model import Lasso, LinearRegression, Ridge
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from sklearn.model_selection import KFold, cross_val_predict, cross_val_score, train_test_split
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

# --------------------------------------------------------------------------
# 1. The dataset -- the only bundled regression set that needs no download
# --------------------------------------------------------------------------


def load_dataset():
    """The Wisconsin/diabetes progression set, in raw measurement units.

    ``load_diabetes`` is the only regression dataset scikit-learn bundles
    offline; ``fetch_california_housing`` downloads and is not used here.
    ``scaled=False`` keeps the ten features in their original units (age in
    years, bmi, average blood pressure, six serum measures) so a coefficient
    would still mean something if this project stopped to interpret one --
    Day 148's point, carried forward. The target itself is a composite
    disease-progression score with no physical unit; it is not measured in
    anything, and this project never pretends otherwise.
    """
    d = load_diabetes(scaled=False)
    return d.data, d.target, list(d.feature_names)


# --------------------------------------------------------------------------
# 2. The frame and the baseline -- before any model
# --------------------------------------------------------------------------


def baseline_metrics(x_train, y_train, x_test, y_test):
    """The mean-predictor baseline: RMSE and R^2, computed before any model.

    Day 141's rule, restated for regression: a score is not evidence until
    you know what it beats. Predicting the training mean for every row is
    the simplest possible non-model.
    """
    dummy = DummyRegressor(strategy="mean").fit(x_train, y_train)
    pred = dummy.predict(x_test)
    rmse = float(np.sqrt(mean_squared_error(y_test, pred)))
    r2 = float(r2_score(y_test, pred))
    return round(rmse, 4), round(r2, 4)


# --------------------------------------------------------------------------
# 3. The split -- 442 rows is small, so the test set stays small too
# --------------------------------------------------------------------------


def split_once(X, y, seed: int = 0, test_size: float = 0.25):
    """One split. The test half is touched once, later, for scoring.

    A 25 percent test set on 442 rows is about 110 rows -- small enough
    that every interval in this project is wide, and that is reported
    honestly rather than smoothed over.
    """
    return train_test_split(X, y, test_size=test_size, random_state=seed)


# --------------------------------------------------------------------------
# 4. The candidate pipelines -- ridge, lasso and plain OLS, K counted
# --------------------------------------------------------------------------

_ALPHAS = [0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30, 100]


def candidate_configs():
    """23 candidate pipelines: 11 ridge, 11 lasso, 1 plain OLS.

    Every candidate is a scikit-learn ``Pipeline`` with a ``StandardScaler``
    ahead of the estimator, so cross-validation refits the scaler on each
    fold's training rows only -- Day 143's stage-ordering rule, enforced by
    the estimator's own contract. Returns ``(family, hyperparameter,
    make_pipeline)`` where ``make_pipeline`` is a zero-argument callable
    returning a fresh, unfitted ``Pipeline``.
    """
    configs = []
    for a in _ALPHAS:
        configs.append(
            ("ridge", a, lambda a=a: Pipeline([("scale", StandardScaler()), ("clf", Ridge(alpha=a))]))
        )
    for a in _ALPHAS:
        configs.append(
            (
                "lasso",
                a,
                lambda a=a: Pipeline(
                    [("scale", StandardScaler()), ("clf", Lasso(alpha=a, max_iter=20000))]
                ),
            )
        )
    configs.append(
        ("ols", 0.0, lambda: Pipeline([("scale", StandardScaler()), ("clf", LinearRegression())]))
    )
    return configs


def candidate_count() -> int:
    """K, the number of configurations this project actually tries."""
    return len(candidate_configs())


# --------------------------------------------------------------------------
# 5. Cross-validate, then select -- honest spending of the train rows
# --------------------------------------------------------------------------


def cross_validate_configs(x_train, y_train, seed: int = 0, folds: int = 5):
    """5-fold CV RMSE for every candidate, on train rows only.

    Returns rows of ``(family, hyperparameter, cv_rmse, cv_std)``, sorted
    best (lowest RMSE) first. RMSE is chosen as the selection metric because
    it is stated as the metric BEFORE any model is fitted -- Day 152's whole
    argument that the metric is a choice, made here in the units the target
    is actually in (unitless composite-score points, not mg/dL).
    """
    splitter = KFold(n_splits=folds, shuffle=True, random_state=seed)
    rows = []
    for family, param, make in candidate_configs():
        scores = cross_val_score(make(), x_train, y_train, cv=splitter, scoring="neg_root_mean_squared_error")
        rows.append((family, param, round(float(-scores.mean()), 4), round(float(scores.std()), 4)))
    rows.sort(key=lambda r: r[2])
    return rows


def select_best(x_train, y_train, seed: int = 0, folds: int = 5):
    """Fit the winner of the sweep on the full training set.

    Returns ``(family, hyperparameter, cv_rmse, fitted_pipeline)``.
    """
    rows = cross_validate_configs(x_train, y_train, seed=seed, folds=folds)
    winner_family, winner_param, winner_cv, _sd = rows[0]
    for family, param, make_fn in candidate_configs():
        if family == winner_family and param == winner_param:
            fitted = make_fn().fit(x_train, y_train)
            return winner_family, winner_param, winner_cv, fitted
    raise RuntimeError("winning configuration vanished between sweep and refit")


# --------------------------------------------------------------------------
# 6. The test set -- one look, enforced mechanically
# --------------------------------------------------------------------------


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


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

    Day 144's discipline, reused unchanged: the counter does not advance on
    a refused attempt, so a caller that never succeeds cannot drain the
    budget by retrying. ``evaluate`` returns ``(rmse, r2, mae)`` -- the one
    look this project spends. Reading predictions afterward for residual
    diagnostics is inspection, not a second selection, exactly as Day 147's
    confusion matrix read predictions after its own one evaluation.
    """

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

    def evaluate(self, model):
        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
        pred = model.predict(self._X)
        rmse = float(np.sqrt(mean_squared_error(self._y, pred)))
        r2 = float(r2_score(self._y, pred))
        mae = float(mean_absolute_error(self._y, pred))
        return round(rmse, 4), round(r2, 4), round(mae, 4)


# --------------------------------------------------------------------------
# 7. The verdict -- a bootstrap interval around the margin
# --------------------------------------------------------------------------


def margin_bootstrap_interval(y_test, pred_baseline, pred_model, n_boot: int = 2000, seed: int = 0):
    """A 95 percent interval around the margin: baseline RMSE minus model RMSE.

    Resamples the test rows with replacement ``n_boot`` times and recomputes
    both RMSEs on each resample, so the interval reflects how much the
    margin itself could plausibly move at this test-set size -- the same
    question Day 144's accuracy interval asked, answered here for RMSE
    without an analytic formula, because RMSE's sampling distribution has
    no equally simple closed form.
    """
    rng = np.random.default_rng(seed)
    y_test = np.asarray(y_test)
    pred_baseline = np.asarray(pred_baseline)
    pred_model = np.asarray(pred_model)
    n = len(y_test)
    margins = np.empty(n_boot)
    for i in range(n_boot):
        idx = rng.integers(0, n, size=n)
        rmse_base = np.sqrt(mean_squared_error(y_test[idx], pred_baseline[idx]))
        rmse_model = np.sqrt(mean_squared_error(y_test[idx], pred_model[idx]))
        margins[i] = rmse_base - rmse_model
    lower = round(float(np.percentile(margins, 2.5)), 4)
    upper = round(float(np.percentile(margins, 97.5)), 4)
    return lower, upper


def margin_distinguishable(lower: float, upper: float) -> bool:
    """Whether the bootstrap interval around the margin excludes zero.

    If the interval spans zero, the honest verdict is "cannot distinguish
    this model from the baseline at this test-set size" -- Day 144's
    cautionary case, now checked for a real margin instead of assumed away.
    """
    return lower > 0.0


# --------------------------------------------------------------------------
# 8. Residual diagnostics -- the centrepiece this day owns
# --------------------------------------------------------------------------


def _normal_ppf(p):
    """The inverse standard-normal CDF, by Acklam's rational approximation.

    Built from scratch because scipy is not one of this lab's three pinned
    dependencies -- Day 153's habit, applied again: when a tool is not
    available, build the minimal piece you need rather than reach past the
    pins. Accurate to about 1.15e-9 across the open interval (0, 1); good
    enough for a Q-Q correlation on 111 points. Public-domain algorithm,
    reference in sources.yml.
    """
    a = [-3.969683028665376e01, 2.209460984245205e02, -2.759285104469687e02,
         1.383577518672690e02, -3.066479806614716e01, 2.506628277459239e00]
    b = [-5.447609879822406e01, 1.615858368580409e02, -1.556989798598866e02,
         6.680131188771972e01, -1.328068155288572e01]
    c = [-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e00,
         -2.549732539343734e00, 4.374664141464968e00, 2.938163982698783e00]
    d = [7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e00,
         3.754408661907416e00]
    p_low = 0.02425
    p = np.asarray(p, dtype=float)
    out = np.empty_like(p)
    low = p < p_low
    high = p > (1 - p_low)
    mid = ~(low | high)
    ql = np.sqrt(-2 * np.log(p[low]))
    out[low] = (((((c[0] * ql + c[1]) * ql + c[2]) * ql + c[3]) * ql + c[4]) * ql + c[5]) / (
        (((d[0] * ql + d[1]) * ql + d[2]) * ql + d[3]) * ql + 1
    )
    q = p[mid] - 0.5
    r = q * q
    out[mid] = (((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * q / (
        (((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1)
    )
    qh = np.sqrt(-2 * np.log(1 - p[high]))
    out[high] = -(((((c[0] * qh + c[1]) * qh + c[2]) * qh + c[3]) * qh + c[4]) * qh + c[5]) / (
        (((d[0] * qh + d[1]) * qh + d[2]) * qh + d[3]) * qh + 1
    )
    return out


def residual_summary(y_test, pred_test):
    """Mean and standard deviation of the test residuals, y minus prediction."""
    resid = np.asarray(y_test, dtype=float) - np.asarray(pred_test, dtype=float)
    return round(float(resid.mean()), 4), round(float(resid.std()), 4)


def heteroscedasticity_signal(pred_test, y_test):
    """Correlation between the fitted value and the absolute residual.

    A value near zero means the spread of errors does not grow with the
    predicted value -- the "fanning out" pattern a residual-vs-fitted plot
    is built to catch. This is the from-scratch numeric version of reading
    that plot.
    """
    pred_test = np.asarray(pred_test, dtype=float)
    resid = np.asarray(y_test, dtype=float) - pred_test
    corr = np.corrcoef(pred_test, np.abs(resid))[0, 1]
    return round(float(corr), 4)


def curvature_signal(pred_test, y_test):
    """Correlation between the squared fitted value and the signed residual.

    A value far from zero suggests the model is missing a systematic curve
    -- residuals that trend up then down (or the reverse) as the fitted
    value rises, rather than scattering evenly around zero.
    """
    pred_test = np.asarray(pred_test, dtype=float)
    resid = np.asarray(y_test, dtype=float) - pred_test
    corr = np.corrcoef(pred_test**2, resid)[0, 1]
    return round(float(corr), 4)


def normal_probability_correlation(y_test, pred_test):
    """A from-scratch normal-probability (Q-Q) check, as one correlation.

    Standardises the residuals, sorts them, and correlates them against the
    theoretical normal quantiles a perfectly Gaussian set of residuals
    would produce. 1.0 is a perfectly straight Q-Q line; values noticeably
    below 1.0 flag departures from normality that a plotted Q-Q line would
    show as curvature at the tails.
    """
    resid = np.asarray(y_test, dtype=float) - np.asarray(pred_test, dtype=float)
    n = len(resid)
    std_resid = (resid - resid.mean()) / resid.std()
    sorted_resid = np.sort(std_resid)
    probs = (np.arange(1, n + 1) - 0.5) / n
    theoretical = _normal_ppf(probs)
    corr = np.corrcoef(theoretical, sorted_resid)[0, 1]
    return round(float(corr), 4)


def largest_residuals(y_test, pred_test, n: int = 5):
    """The n largest-magnitude residuals, inspected individually.

    Returns rows of ``(test_row_index, true_value, predicted_value,
    residual)``, sorted by |residual| descending. A confusion matrix reads
    the specific mistakes a classifier makes; this is that discipline's
    regression counterpart.
    """
    y_test = np.asarray(y_test, dtype=float)
    pred_test = np.asarray(pred_test, dtype=float)
    resid = y_test - pred_test
    order = np.argsort(-np.abs(resid))[:n]
    return [
        (int(i), round(float(y_test[i]), 4), round(float(pred_test[i]), 4), round(float(resid[i]), 4))
        for i in order
    ]


# --------------------------------------------------------------------------
# 9. Is the model worse for high-value targets? Measure it.
# --------------------------------------------------------------------------


def error_by_target_level(y_test, pred_test):
    """RMSE on the below-median half of test targets against the above-median half.

    Returns ``(rmse_low, rmse_high, ratio)`` where ``ratio`` is
    ``rmse_high / rmse_low``. On a disease-progression score this is a
    fairness-relevant question, not only a statistical one: are errors
    worse for patients whose true progression is more severe?
    """
    y_test = np.asarray(y_test, dtype=float)
    pred_test = np.asarray(pred_test, dtype=float)
    median_y = float(np.median(y_test))
    low_mask = y_test <= median_y
    high_mask = ~low_mask
    rmse_low = float(np.sqrt(mean_squared_error(y_test[low_mask], pred_test[low_mask])))
    rmse_high = float(np.sqrt(mean_squared_error(y_test[high_mask], pred_test[high_mask])))
    return round(rmse_low, 4), round(rmse_high, 4), round(rmse_high / rmse_low, 4)


# --------------------------------------------------------------------------
# 10. The leaky version -- selecting by peeking at the test set
# --------------------------------------------------------------------------


def leaky_selection_test_rmse(x_train, y_train, x_test, y_test):
    """Select the winner by fitting every candidate and scoring it on TEST.

    Returns the winning (lowest) test RMSE -- K looks disguised as one,
    Day 147's mistake, reconstructed for regression. Lower RMSE is better,
    so the leaky search can only match or beat the honestly selected
    model's own test RMSE, never lose to it.
    """
    best_rmse = None
    for _family, _param, make in candidate_configs():
        pipe = make().fit(x_train, y_train)
        pred = pipe.predict(x_test)
        rmse = float(np.sqrt(mean_squared_error(y_test, pred)))
        if best_rmse is None or rmse < best_rmse:
            best_rmse = rmse
    return round(best_rmse, 4)


def leaky_vs_honest_over_seeds(X, y, seeds=range(20), folds: int = 5):
    """The gap between peeking at the test set and looking at it once.

    Returns rows of ``(seed, honest_rmse, leaky_rmse, gap)`` where
    ``gap = honest_rmse - leaky_rmse``. Because lower RMSE is better, a
    positive gap means the leak reported a lower (better-looking) error
    than the honest evaluation -- the leak can only help the reported
    number, never hurt it.
    """
    rows = []
    for seed in seeds:
        x_train, x_test, y_train, y_test = split_once(X, y, seed=seed)
        _family, _param, _cv, fitted = select_best(x_train, y_train, seed=seed, folds=folds)
        honest_rmse = float(np.sqrt(mean_squared_error(y_test, fitted.predict(x_test))))
        leaky_rmse = leaky_selection_test_rmse(x_train, y_train, x_test, y_test)
        rows.append((seed, round(honest_rmse, 4), leaky_rmse, round(honest_rmse - leaky_rmse, 4)))
    return rows


# --------------------------------------------------------------------------
# 11. Prediction intervals -- not just a point, and measured coverage
# --------------------------------------------------------------------------


def prediction_interval_coverage(x_train, y_train, x_test, y_test, fitted, seed: int = 0, folds: int = 5):
    """A constant-width 95 percent prediction interval, and its realised coverage.

    The half-width comes from the standard deviation of out-of-fold
    residuals on the TRAINING rows only (``cross_val_predict``), never from
    the test residuals themselves -- using the test residuals to size the
    test interval would be circular. Returns ``(half_width, coverage)``
    where ``coverage`` is the fraction of test targets that actually fall
    inside ``prediction +/- half_width``, measured against the nominal 0.95.
    """
    splitter = KFold(n_splits=folds, shuffle=True, random_state=seed)
    oof_pred = cross_val_predict(fitted, x_train, y_train, cv=splitter)
    oof_resid = np.asarray(y_train, dtype=float) - oof_pred
    half_width = round(float(1.96 * oof_resid.std()), 4)
    pred_test = fitted.predict(x_test)
    y_test = np.asarray(y_test, dtype=float)
    lower = pred_test - half_width
    upper = pred_test + half_width
    coverage = round(float(np.mean((y_test >= lower) & (y_test <= upper))), 4)
    return half_width, coverage

"""Fourteen exercises: one regression project, run properly, once.

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.
`regression_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 inspect

import numpy as np  # noqa: F401  (you will need it)
import pytest
from sklearn.datasets import fetch_california_housing, load_diabetes  # noqa: F401
from sklearn.dummy import DummyRegressor  # noqa: F401

import regression_lib as r  # noqa: F401  (you will need it)


@pytest.fixture(scope="module")
def dataset():
    return r.load_dataset()


@pytest.fixture(scope="module")
def split(dataset):
    X, y, _names = dataset
    return r.split_once(X, y, seed=0)


def test_01_the_dataset_is_the_only_bundled_regression_set(dataset):
    pytest.skip(
        "Unpack (X, y, names) from the dataset fixture. Assert X.shape == "
        "(442, 10), names == ['age', 'sex', 'bmi', 'bp', 's1', 's2', 's3', "
        "'s4', 's5', 's6'], y.min() rounds to 25.0, y.max() rounds to 346.0, "
        "and y.mean() rounds to 152.1335. Then assert "
        "inspect.signature(fetch_california_housing).parameters"
        "['download_if_missing'].default is True -- the reason that dataset "
        "is forbidden here and load_diabetes is used instead."
    )


def test_01b_raw_units_are_not_the_scikit_learn_default(dataset):
    pytest.skip(
        "Unpack (X, y, names) from the dataset fixture. Find age's column "
        "index and assert X[:, age_index].min() == 19.0 and .max() == 79.0 "
        "-- real years. Then load load_diabetes(scaled=True) (scikit-learn's "
        "own default) and assert its age column's min is between -0.2 and 0, "
        "and its max is between 0 and 0.2 -- tiny centred floats. This is "
        "why this lab asks for scaled=False."
    )


def test_02_the_baseline_before_any_model(split):
    pytest.skip(
        "Unpack (x_train, x_test, y_train, y_test) from the split fixture. "
        "Call r.baseline_metrics and assert the RMSE rounds to 70.4637 and "
        "R2 rounds to -0.0001. Every model in this exercise has to beat "
        "70.4637 RMSE to be worth building at all."
    )


def test_03_the_split_holds_the_test_rows_back(dataset):
    pytest.skip(
        "Unpack (X, y, _names) from the dataset fixture and call "
        "r.split_once(X, y, seed=0). Assert x_train.shape == (331, 10) and "
        "x_test.shape == (111, 10). 111 test rows on 442 total is small -- "
        "every interval computed later in this exercise is wide because of "
        "it, and that is the point, not a bug."
    )


def test_04_the_sweep_counts_twenty_three_candidate_pipelines():
    pytest.skip(
        "Assert r.candidate_count() == 23. Then unpack the (family, param, "
        "make) triples from r.candidate_configs() and assert there are 11 "
        "'ridge', 11 'lasso' and 1 'ols' entries. K is the number nobody "
        "remembers -- count it before doing anything else with it."
    )


def test_05_cross_validation_selects_the_winner_on_train_rows_only(split):
    pytest.skip(
        "Unpack (x_train, _x_test, y_train, _y_test) from the split fixture "
        "and call r.select_best(x_train, y_train, seed=0). Assert the "
        "returned (family, param) equals ('lasso', 1) and cv_rmse equals "
        "53.8958. The winner was chosen on cross-validated train rows -- "
        "the test rows have not been touched yet."
    )


def test_06_the_gate_permits_exactly_one_test_evaluation(split):
    pytest.skip(
        "Fit the winner from r.select_best on the train rows, wrap "
        "(x_test, y_test) in r.GatedTestSet, and assert the first "
        "evaluation returns (56.5566, 0.3557, 45.2846) for (rmse, r2, mae) "
        "and the counter becomes 1. Then assert a second evaluation raises "
        "r.TestSetTouchedTwice mentioning 'validation score', and that the "
        "counter did NOT advance on the refused attempt."
    )


def test_07_the_margin_has_a_bootstrap_interval(split):
    pytest.skip(
        "Fit a DummyRegressor(strategy='mean') baseline and the winner from "
        "r.select_best, both on x_train/y_train. Predict both on x_test. "
        "Call r.margin_bootstrap_interval(y_test, pred_baseline, pred_model, "
        "seed=0) and assert it returns (5.5852, 22.3324). Assert the point "
        "margin -- baseline RMSE minus model RMSE -- rounds to 13.9071, and "
        "that r.margin_distinguishable(lower, upper) is True: the margin "
        "clears the interval's lower bound."
    )


def test_08_the_residual_vs_fitted_diagnostic(split):
    pytest.skip(
        "Fit the winner and predict on x_test. Call r.residual_summary and "
        "assert it returns (-3.6262, 56.4402). Call "
        "r.heteroscedasticity_signal(pred, y_test) and assert it rounds to "
        "0.2386. Call r.curvature_signal(pred, y_test) and assert it rounds "
        "to -0.1278. Neither signal is dramatic here -- report both anyway."
    )


def test_08b_the_normal_probability_check_and_the_largest_residuals(split):
    pytest.skip(
        "Fit the winner and predict on x_test. Call "
        "r.normal_probability_correlation(y_test, pred) and assert it "
        "rounds to 0.9901. Call r.largest_residuals(y_test, pred, n=5) and "
        "assert the first two rows equal (60, 52.0, 209.3314, -157.3314) "
        "and (65, 302.0, 153.8865, 148.1135). 0.9901 is close to a straight "
        "Q-Q line -- these residuals are close to normal even on real data."
    )


def test_09_error_by_target_level(split):
    pytest.skip(
        "Fit the winner and predict on x_test. Call "
        "r.error_by_target_level(y_test, pred) and assert it returns "
        "(55.2464, 57.8601, 1.0473). A ratio of 1.0473 means the model is "
        "only 4.73 percent worse on the more-severe half of test targets -- "
        "measure it before assuming a disease-progression model is fair or "
        "unfair across severity."
    )


def test_10_the_leaky_version_selects_by_peeking_at_the_test_set(split):
    pytest.skip(
        "Fit the honest winner and compute its test RMSE (round to 4 "
        "places). Call r.leaky_selection_test_rmse(x_train, y_train, "
        "x_test, y_test), which fits every one of the 23 candidates and "
        "lets the test set itself pick the winner by lowest RMSE. Assert "
        "the honest RMSE is 56.5566, the leaky RMSE is 55.5212, and "
        "leaky_rmse <= honest_rmse -- lower RMSE is better, so the leak can "
        "only match or beat the honest score, never lose to it."
    )


def test_10b_the_leaky_gap_over_twenty_seeds(dataset):
    pytest.skip(
        "Call r.leaky_vs_honest_over_seeds(X, y, seeds=range(20)). Assert "
        "20 rows come back. Compute the mean, sd, min and max of the gap "
        "column and assert they round to 0.5279, 0.3686, 0.011 and 1.1451. "
        "Then assert every single gap is non-negative -- across 20 "
        "independent seeds, selecting by peeking at the test set never once "
        "reported a worse (higher) RMSE than selecting honestly."
    )


def test_11_prediction_interval_coverage(split):
    pytest.skip(
        "Fit the winner. Call r.prediction_interval_coverage(x_train, "
        "y_train, x_test, y_test, fitted, seed=0) and assert it returns "
        "(105.8797, 0.9459). The half-width comes from TRAINING residuals "
        "only, never from the test residuals themselves -- sizing an "
        "interval from the data you are about to check it against would be "
        "circular. 0.9459 measured against a 0.95 nominal, on 111 test rows."
    )

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

These five 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 regression_lib as r


def test_the_dataset_loads_offline_and_matches_its_shape():
    X, y, names = r.load_dataset()
    assert X.shape == (442, 10)
    assert y.shape == (442,)
    assert names == ["age", "sex", "bmi", "bp", "s1", "s2", "s3", "s4", "s5", "s6"]
    assert float(y.min()) == 25.0
    assert float(y.max()) == 346.0


def test_candidate_configs_really_are_twenty_three_distinct_pipelines():
    configs = r.candidate_configs()
    assert len(configs) == 23
    seen = set()
    for family, param, make in configs:
        pipe = make()
        assert hasattr(pipe, "fit") and hasattr(pipe, "predict")
        seen.add((family, param))
    # No two configs share a (family, hyperparameter) pair.
    assert len(seen) == 23


def test_the_gated_test_set_counts_and_refuses():
    class AlwaysMean:
        def predict(self, X):
            return np.full(len(X), 150.0)

    y = np.array([100.0, 200.0, 150.0, 150.0])
    gate = r.GatedTestSet(np.zeros((4, 2)), y)
    rmse, r2, mae = gate.evaluate(AlwaysMean())
    assert rmse > 0
    with pytest.raises(r.TestSetTouchedTwice):
        gate.evaluate(AlwaysMean())
    # A fresh gate is a fresh budget; the class holds no global state.
    fresh_rmse, _r2, _mae = r.GatedTestSet(np.zeros((4, 2)), y).evaluate(AlwaysMean())
    assert fresh_rmse == rmse


def test_the_normal_ppf_matches_known_reference_points():
    # Well-known standard-normal quantiles, to a few decimals.
    assert abs(r._normal_ppf(np.array([0.5]))[0] - 0.0) < 1e-6
    assert abs(r._normal_ppf(np.array([0.975]))[0] - 1.959964) < 1e-4
    assert abs(r._normal_ppf(np.array([0.025]))[0] - (-1.959964)) < 1e-4


def test_largest_residuals_returns_them_sorted_by_magnitude():
    y_test = np.array([10.0, 20.0, 30.0, 40.0])
    pred_test = np.array([10.0, 25.0, 10.0, 41.0])
    rows = r.largest_residuals(y_test, pred_test, n=2)
    assert len(rows) == 2
    # Row 2 (residual +20) is the largest in magnitude; row 3 (residual -1) is smallest.
    assert rows[0][0] == 2
    assert rows[0][3] == 20.0

#!/usr/bin/env bash
# Day 154 lab harness: "A Complete Regression Project"
#
# 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.dummy import DummyRegressor

import regression_lib as r

errors = []


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


# 1. The dataset
X, y, names = r.load_dataset()
expect("shape", X.shape, (442, 10))
expect("y min", round(float(y.min()), 4), 25.0)
expect("y max", round(float(y.max()), 4), 346.0)
expect("y mean", round(float(y.mean()), 4), 152.1335)

# 2. The split and the baseline
x_train, x_test, y_train, y_test = r.split_once(X, y, seed=0)
expect("train shape", x_train.shape, (331, 10))
expect("test shape", x_test.shape, (111, 10))
base_rmse, base_r2 = r.baseline_metrics(x_train, y_train, x_test, y_test)
expect("baseline rmse", base_rmse, 70.4637)
expect("baseline r2", base_r2, -0.0001)

# 3. The sweep and selection
expect("K", r.candidate_count(), 23)
family, param, cv_rmse, fitted = r.select_best(x_train, y_train, seed=0)
expect("winner", (family, param), ("lasso", 1))
expect("cv_rmse", cv_rmse, 53.8958)

# 4. The gate: exactly one evaluation
gate = r.GatedTestSet(x_test, y_test)
test_rmse, test_r2, test_mae = gate.evaluate(fitted)
expect("test_rmse", test_rmse, 56.5566)
expect("test_r2", test_r2, 0.3557)
expect("test_mae", test_mae, 45.2846)
expect("evaluations after first look", gate.evaluations, 1)
try:
    gate.evaluate(fitted)
    errors.append("the gate permitted a second evaluation, which it must not")
except r.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")

# 5. The margin, with a bootstrap interval
baseline_model = DummyRegressor(strategy="mean").fit(x_train, y_train)
pred_baseline = baseline_model.predict(x_test)
pred_model = fitted.predict(x_test)
lower, upper = r.margin_bootstrap_interval(y_test, pred_baseline, pred_model, seed=0)
expect("margin interval", (lower, upper), (5.5852, 22.3324))
margin = round(base_rmse - test_rmse, 4)
expect("margin", margin, 13.9071)
expect("distinguishable", r.margin_distinguishable(lower, upper), True)

# 6. Residual diagnostics
resid_mean, resid_std = r.residual_summary(y_test, pred_model)
expect("resid_mean", resid_mean, -3.6262)
expect("resid_std", resid_std, 56.4402)
expect("heteroscedasticity", r.heteroscedasticity_signal(pred_model, y_test), 0.2386)
expect("curvature", r.curvature_signal(pred_model, y_test), -0.1278)
expect("qq correlation", r.normal_probability_correlation(y_test, pred_model), 0.9901)
top5 = r.largest_residuals(y_test, pred_model, n=5)
expect("largest residual row", top5[0], (60, 52.0, 209.3314, -157.3314))

# 7. Error by target level
rmse_low, rmse_high, ratio = r.error_by_target_level(y_test, pred_model)
expect("error by level", (rmse_low, rmse_high, ratio), (55.2464, 57.8601, 1.0473))

# 8. The leaky version
leaky_rmse = r.leaky_selection_test_rmse(x_train, y_train, x_test, y_test)
expect("leaky rmse", leaky_rmse, 55.5212)
if leaky_rmse > test_rmse:
    errors.append(f"leaky rmse {leaky_rmse} was worse (higher) than the honest rmse {test_rmse}")

# 9. Prediction interval coverage
half_width, coverage = r.prediction_interval_coverage(x_train, y_train, x_test, y_test, fitted, seed=0)
expect("half width", half_width, 105.8797)
expect("coverage", coverage, 0.9459)

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-11 reproduced directly against regression_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 "^19 passed"; then
  ok "pytest examples -q -> 19 passed"
else
  fail "pytest examples -q did not report 19 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 "5 passed, 14 skipped"; then
  ok "pytest starter -q -> 5 passed, 14 skipped (the machinery checks pass; the fourteen exercises are stubs)"
else
  fail "pytest starter -q did not report 5 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. The test set is evaluated EXACTLY ONCE in the reference run"
GATE_CHECK=$("$PYTHON" - <<'PYEOF'
import sys
sys.path.insert(0, "examples")
import regression_lib as r

X, y, _names = r.load_dataset()
x_train, x_test, y_train, y_test = r.split_once(X, y, seed=0)
_family, _param, _cv, fitted = r.select_best(x_train, y_train, seed=0)

gate = r.GatedTestSet(x_test, y_test)
assert gate.evaluations == 0, "a fresh gate must start at zero evaluations"
gate.evaluate(fitted)
assert gate.evaluations == 1, "one evaluation must advance the counter to exactly one"
refused = 0
for _ in range(5):
    try:
        gate.evaluate(fitted)
    except r.TestSetTouchedTwice:
        refused += 1
assert refused == 5, "every repeated attempt after the first must be refused"
assert gate.evaluations == 1, "repeated refused attempts must not advance the counter"
print("test set touched exactly once, five further attempts refused, counter never moved")
PYEOF
)
if echo "$GATE_CHECK" | grep -q "touched exactly once"; then
  ok "GatedTestSet enforces exactly one evaluation mechanically, not by convention"
else
  fail "the one-evaluation guarantee did not hold"
  echo "$GATE_CHECK" | sed 's/^/    /'
fi

echo ""
echo "8. Proof the harness can fail"
SCRATCH=$(mktemp -d "${TMPDIR:-/tmp}/d154-scratch.XXXXXX")
cp examples/*.py "$SCRATCH"/
SCRATCH_OUT=$("$PYTEST" "$SCRATCH" -q 2>&1)
if echo "$SCRATCH_OUT" | tail -1 | grep -qE "^19 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_regression_claims.py" <<'PYEOF'
import sys
path = sys.argv[1]
text = open(path).read()
needle = "assert rows[0] == (60, 52.0, 209.3314, -157.3314)"
replacement = "assert rows[0] == (0, 0.0, 0.0, 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_08b_the_normal_probability_check_and_the_largest_residuals"; then
  ok "breaking exercise 8b'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 "9. The leaky-gap direction holds beyond the quoted seed range"
DIRECTION=$("$PYTHON" - <<'PYEOF'
import sys
sys.path.insert(0, "examples")
import numpy as np
import regression_lib as r

X, y, _names = r.load_dataset()
problems = []

# The leak never produces a worse (higher) RMSE, at seeds this lab does not quote.
rows = r.leaky_vs_honest_over_seeds(X, y, seeds=range(20, 25))
gaps = [g for _s, _h, _l, g in rows]
if not all(g >= 0 for g in gaps):
    problems.append(f"a leaky gap went negative at an unquoted seed: {gaps}")

# Selecting is confined to train rows: refitting the winner never needs test.
x_train, x_test, y_train, y_test = r.split_once(X, y, seed=41)
_family, _param, cv_rmse, fitted = r.select_best(x_train, y_train, seed=41)
if not (30 < cv_rmse < 90):
    problems.append(f"cv_rmse at an unquoted seed was out of range: {cv_rmse}")

if problems:
    for p in problems:
        print("ERROR:", p)
else:
    print("every direction held")
PYEOF
)
if [ "$DIRECTION" = "every direction held" ]; then
  ok "the leaky-gap direction and the selection mechanics hold at seeds this lab does not quote"
else
  fail "a direction failed beyond the quoted seeds"
  echo "$DIRECTION" | sed 's/^/    /'
fi

echo ""
echo "10. 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 beyond the bundled dataset"
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