Hands-on lab — Day 150: Multiple and Polynomial Regression

"""Many predictors, measured: what changes when there is more than one.

Day 148 covered the line through one predictor. Day 149 covered why the
loss is squared error. This module measures what is new once a second
predictor joins the first: a coefficient's meaning becomes conditional on
"holding the others constant", and when predictors are correlated with
each other -- not just with the target -- that condition gets expensive.

The centrepiece is a duplicated predictor. Two columns that carry the same
information cannot be told apart by the normal equations, so the fit
smears one true effect across both coefficients in whatever proportion
happens to minimise squared error that day -- while the *sum* of the two
coefficients, and every prediction the model makes, barely moves at all.

Everything here is deterministic given a seed, and every dataset is either
the bundled `sklearn.datasets.load_diabetes(scaled=False)` -- ten raw-unit
clinical predictors, real measurements, no download -- or built from it.
"""

from __future__ import annotations

import numpy as np

from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures, StandardScaler


# --------------------------------------------------------------------------
# 0. The dataset itself
# --------------------------------------------------------------------------


def load_raw_diabetes():
    """The ten raw-unit predictors and the target, with feature names.

    `scaled=False` matters: it returns age in years, sex coded 1/2, bmi,
    average blood pressure, and six serum measurements s1-s6 in their
    original units, which is what makes a coefficient's magnitude mean
    anything. The default `scaled=True` mean-centres and unit-norm scales
    every column before you ever see it -- exercise 8 measures exactly
    what that substitution does and does not change.
    """
    bunch = load_diabetes(scaled=False)
    return bunch.data, bunch.target, list(bunch.feature_names)


def fit(X, y):
    """A plain ordinary-least-squares fit -- the one tool this lesson uses."""
    return LinearRegression().fit(X, y)


# --------------------------------------------------------------------------
# 1. Correlation and variance inflation
# --------------------------------------------------------------------------


def correlation(X, names, a: str, b: str) -> float:
    """The Pearson correlation between two named predictor columns."""
    i, j = names.index(a), names.index(b)
    return float(np.corrcoef(X[:, i], X[:, j])[0, 1])


def variance_inflation_factors(X, names) -> dict:
    """VIF for every column: regress it on the rest, VIF = 1 / (1 - R2).

    A VIF of 1 means a predictor is unrelated to the others. Above 5 or 10
    is the usual rule of thumb for "correlated enough to worry about" --
    this function does not apply the rule, it only computes the number the
    rule is applied to.
    """
    result = {}
    for i, name in enumerate(names):
        others = np.delete(X, i, axis=1)
        target = X[:, i]
        r2 = LinearRegression().fit(others, target).score(others, target)
        result[name] = float("inf") if r2 >= 1.0 else round(1.0 / (1.0 - r2), 4)
    return result


# --------------------------------------------------------------------------
# 2. The centrepiece: an exact duplicate column
# --------------------------------------------------------------------------


def duplicate_column_exact(X, y, col_index: int):
    """Append an exact copy of one column and refit.

    Returns ``(original_coef, dup_coef_a, dup_coef_b, max_abs_pred_diff,
    r2_original, r2_dup)``. The two duplicate coefficients need not equal
    the original one individually -- only their sum does, because the
    normal equations only ever "see" the combined effect of two identical
    columns, and there is no unique way to split it.
    """
    original = fit(X, y)
    X_dup = np.hstack([X, X[:, [col_index]]])
    dup = fit(X_dup, y)
    pred_original = original.predict(X)
    pred_dup = dup.predict(X_dup)
    return (
        float(original.coef_[col_index]),
        float(dup.coef_[col_index]),
        float(dup.coef_[-1]),
        float(np.max(np.abs(pred_dup - pred_original))),
        float(original.score(X, y)),
        float(dup.score(X_dup, y)),
    )


def duplicate_column_noisy(X, y, col_index: int, noise_scale: float, seed: int):
    """Append a near-duplicate: the same column plus a little noise.

    Breaking the exact tie lets ordinary least squares pick a *unique*
    split of the combined effect again -- but which split it picks depends
    on which way the noise happened to fall, which is the whole point.
    Returns ``(coef_a, coef_b, sum_coefs, max_abs_pred_diff, r2)``.
    """
    original = fit(X, y)
    pred_original = original.predict(X)
    rng = np.random.default_rng(seed)
    noise = rng.normal(scale=noise_scale, size=X.shape[0])
    near_dup = X[:, col_index] + noise
    X_noisy = np.hstack([X, near_dup.reshape(-1, 1)])
    model = fit(X_noisy, y)
    pred = model.predict(X_noisy)
    return (
        float(model.coef_[col_index]),
        float(model.coef_[-1]),
        float(model.coef_[col_index] + model.coef_[-1]),
        float(np.max(np.abs(pred - pred_original))),
        float(model.score(X_noisy, y)),
    )


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


def duplicate_noisy_spread(X, y, col_index: int, noise_scale: float, seeds) -> dict:
    """Across many noise draws: how the two coefficients wander, and how
    their sum, the predictions and R2 do not.

    Returns a dict with ``coef_a``, ``coef_b`` and ``sum`` spreads (each
    from :func:`spread`), plus ``max_pred_diff`` (the largest single
    prediction movement seen across every seed) and an ``r2`` spread.
    """
    coef_a, coef_b, coef_sum, max_diffs, r2s = [], [], [], [], []
    for seed in seeds:
        a, b, total, max_diff, r2 = duplicate_column_noisy(X, y, col_index, noise_scale, seed)
        coef_a.append(a)
        coef_b.append(b)
        coef_sum.append(total)
        max_diffs.append(max_diff)
        r2s.append(r2)
    return {
        "coef_a": spread(coef_a),
        "coef_b": spread(coef_b),
        "sum": spread(coef_sum),
        "max_pred_diff_overall": round(max(max_diffs), 4),
        "r2": spread(r2s),
    }


# --------------------------------------------------------------------------
# 3. Bootstrap coefficient instability, across all ten predictors
# --------------------------------------------------------------------------


def bootstrap_coefficient_spread(X, y, names, reps: int = 500, seed: int = 0) -> dict:
    """Resample the rows with replacement, refit, and record each coefficient.

    Returns ``{name: {"mean": ..., "sd": ..., "cv": ...}}`` where ``cv`` is
    the coefficient of variation, ``abs(sd / mean)`` -- a scale-free way to
    compare how much a coefficient wanders relative to its own size.
    """
    rng = np.random.default_rng(seed)
    n = X.shape[0]
    collected = {name: [] for name in names}
    for _ in range(reps):
        idx = rng.integers(0, n, size=n)
        model = fit(X[idx], y[idx])
        for i, name in enumerate(names):
            collected[name].append(model.coef_[i])
    result = {}
    for name in names:
        arr = np.asarray(collected[name])
        mean = float(arr.mean())
        sd = float(arr.std())
        result[name] = {
            "mean": round(mean, 4),
            "sd": round(sd, 4),
            "cv": round(abs(sd / mean), 4) if mean != 0 else None,
        }
    return result


# --------------------------------------------------------------------------
# 4. What "holding the others constant" changes: simple vs. multiple
# --------------------------------------------------------------------------


def simple_vs_multiple_coefficients(X, y, names) -> dict:
    """Each predictor's coefficient alone, and again inside the full model.

    Returns ``{name: {"simple": ..., "multiple": ..., "sign_flip": bool}}``.
    A sign flip means the predictor's *apparent* relationship with the
    target, taken alone, points the opposite way from its *conditional*
    relationship once the other nine predictors are held constant.
    """
    result = {}
    full = fit(X, y)
    for i, name in enumerate(names):
        alone = fit(X[:, [i]], y)
        simple_coef = float(alone.coef_[0])
        multiple_coef = float(full.coef_[i])
        result[name] = {
            "simple": round(simple_coef, 4),
            "multiple": round(multiple_coef, 4),
            "sign_flip": bool(np.sign(simple_coef) != np.sign(multiple_coef)),
        }
    return result


# --------------------------------------------------------------------------
# 5. A polynomial fit is linear in its parameters
# --------------------------------------------------------------------------


def polynomial_matches_normal_equations(X2, y, degree: int = 2, feature_names=None):
    """PolynomialFeatures + LinearRegression against a direct normal-equations solve.

    Builds the degree-``degree`` design matrix two ways: once through
    scikit-learn's transformer-plus-estimator pipeline, once by hand with
    ``numpy.linalg.lstsq`` on the identical expanded matrix. If a
    polynomial fit really is "linear in its parameters, not in x", the two
    must agree to floating-point precision, because they are solving the
    same linear system.

    Returns ``(feature_names, sklearn_coefs, sklearn_intercept,
    normal_eq_coefs, normal_eq_intercept, max_abs_coef_diff,
    max_abs_intercept_diff)``.
    """
    poly = PolynomialFeatures(degree=degree, include_bias=False)
    X_poly = poly.fit_transform(X2)
    expanded_names = list(poly.get_feature_names_out(feature_names))

    sklearn_model = fit(X_poly, y)

    design = np.hstack([np.ones((X_poly.shape[0], 1)), X_poly])
    beta, *_ = np.linalg.lstsq(design, y, rcond=None)

    coef_diff = float(np.max(np.abs(beta[1:] - sklearn_model.coef_)))
    intercept_diff = float(abs(beta[0] - sklearn_model.intercept_))

    return (
        expanded_names,
        [round(float(c), 6) for c in sklearn_model.coef_],
        round(float(sklearn_model.intercept_), 6),
        [round(float(b), 6) for b in beta[1:]],
        round(float(beta[0]), 6),
        coef_diff,
        intercept_diff,
    )


def interaction_term_effect(X2, y):
    """Compare a degree-2 fit with and without its interaction term.

    ``bmi^2`` and ``bp^2`` describe how each predictor curves on its own.
    ``bmi bp`` describes something neither can: whether the *effect* of
    one depends on the level of the other. Drop it and refit on the
    remaining four columns; the gap in R2 is what the interaction term was
    buying.

    Returns ``(r2_with_interaction, r2_without_interaction,
    interaction_coefficient)``.
    """
    poly = PolynomialFeatures(degree=2, include_bias=False)
    X_poly = poly.fit_transform(X2)
    feature_names = list(poly.get_feature_names_out(["a", "b"]))
    interaction_index = feature_names.index("a b")

    with_interaction = fit(X_poly, y)
    X_without = np.delete(X_poly, interaction_index, axis=1)
    without_interaction = fit(X_without, y)

    return (
        float(with_interaction.score(X_poly, y)),
        float(without_interaction.score(X_without, y)),
        float(with_interaction.coef_[interaction_index]),
    )


# --------------------------------------------------------------------------
# 6. R-squared never decreases when you add a predictor -- even noise
# --------------------------------------------------------------------------


def r2_with_added_noise_columns(X, y, noise_counts, seed: int = 42) -> list:
    """R2 of the fit as pure-noise columns are appended, one count at a time.

    Returns rows of ``(n_noise_columns, r2, delta_from_baseline)``. Every
    added column is `numpy.random.default_rng` noise with no relationship
    to the target whatsoever -- Day 152 owns the fix (adjusted R2); this
    function only measures the problem it fixes.
    """
    base_r2 = fit(X, y).score(X, y)
    rng = np.random.default_rng(seed)
    rows = []
    for n_noise in noise_counts:
        noise_columns = rng.normal(size=(X.shape[0], n_noise))
        X_aug = np.hstack([X, noise_columns])
        r2 = fit(X_aug, y).score(X_aug, y)
        rows.append((n_noise, round(r2, 6), round(r2 - base_r2, 6)))
    return rows


# --------------------------------------------------------------------------
# 7. Scaling changes the coefficients, not the model
# --------------------------------------------------------------------------


def scaling_effect(X, y):
    """Fit on raw units and on standardised units; compare everything.

    Returns ``(raw_coefs, scaled_coefs, raw_r2, scaled_r2,
    max_abs_pred_diff)``. Standardising centres and unit-variance-scales
    every column before fitting, which changes what one unit of a
    predictor means and therefore changes every coefficient's size -- but
    changes nothing about what the model actually predicts.
    """
    raw_model = fit(X, y)
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)
    scaled_model = fit(X_scaled, y)

    pred_raw = raw_model.predict(X)
    pred_scaled = scaled_model.predict(X_scaled)

    return (
        [round(float(c), 4) for c in raw_model.coef_],
        [round(float(c), 4) for c in scaled_model.coef_],
        float(raw_model.score(X, y)),
        float(scaled_model.score(X_scaled, y)),
        float(np.max(np.abs(pred_raw - pred_scaled))),
    )

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

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

import sys
from pathlib import Path

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

import numpy as np  # noqa: E402

import regression_lib as r  # noqa: E402


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


def main() -> None:
    print("Day 150 -- multiple and polynomial regression, measured")
    print("=" * 55)

    X, y, names = r.load_raw_diabetes()

    rule("1. Correlation and variance inflation, ten raw-unit predictors")
    vifs = r.variance_inflation_factors(X, names)
    print("  name    VIF")
    for name in names:
        print(f"  {name:<5}   {vifs[name]:>8.4f}")
    print(f"  correlation(s1, s2) : {r.correlation(X, names, 's1', 's2'):+.4f}")
    print(f"  correlation(s3, s4) : {r.correlation(X, names, 's3', 's4'):+.4f}")
    print(f"  correlation(bmi, bp): {r.correlation(X, names, 'bmi', 'bp'):+.4f}")

    rule("2. The centrepiece: an exact duplicate of s1")
    idx_s1 = names.index("s1")
    original, coef_a, coef_b, max_diff, r2_orig, r2_dup = r.duplicate_column_exact(X, y, idx_s1)
    print(f"  original s1 coefficient        : {original:+.4f}")
    print(f"  duplicate model, coefficient a : {coef_a:+.4f}")
    print(f"  duplicate model, coefficient b : {coef_b:+.4f}")
    print(f"  sum of the two                 : {coef_a + coef_b:+.4f}")
    print(f"  max abs prediction difference  : {max_diff:.2e}")
    print(f"  R2 original / duplicate        : {r2_orig:.4f} / {r2_dup:.4f}")

    rule("3. Breaking the tie with noise")
    noise_scale = 0.01 * float(X[:, idx_s1].std())
    coef_a, coef_b, coef_sum, max_diff, r2 = r.duplicate_column_noisy(X, y, idx_s1, noise_scale, seed=0)
    print(f"  noise scale (1% of s1's std)   : {noise_scale:.4f}")
    print(f"  seed 0, coefficient a          : {coef_a:+.4f}")
    print(f"  seed 0, coefficient b          : {coef_b:+.4f}")
    print(f"  seed 0, sum                    : {coef_sum:+.4f}")
    print(f"  seed 0, R2                     : {r2:.4f}")
    spread10 = r.duplicate_noisy_spread(X, y, idx_s1, noise_scale, range(10))
    print(f"  across seeds 0-9, coefficient a : {spread10['coef_a']}")
    print(f"  across seeds 0-9, coefficient b : {spread10['coef_b']}")
    print(f"  across seeds 0-9, sum           : {spread10['sum']}")
    print(f"  largest single prediction move : {spread10['max_pred_diff_overall']:.4f}")
    print(f"  across seeds 0-9, R2            : {spread10['r2']}")

    rule("4. Bootstrap coefficient instability, all ten predictors")
    boot = r.bootstrap_coefficient_spread(X, y, names, reps=500, seed=0)
    print("  name    mean       sd        cv")
    for name in names:
        row = boot[name]
        cv = "  n/a " if row["cv"] is None else f"{row['cv']:.4f}"
        print(f"  {name:<5}   {row['mean']:>8.4f}   {row['sd']:>7.4f}   {cv}")

    rule("5. Holding the other nine predictors constant flips a sign")
    svm = r.simple_vs_multiple_coefficients(X, y, names)
    print("  name    simple      multiple    sign flip")
    for name in names:
        row = svm[name]
        print(f"  {name:<5}   {row['simple']:>9.4f}   {row['multiple']:>9.4f}   {row['sign_flip']}")

    rule("6. A polynomial fit is linear in its parameters")
    idx_bmi, idx_bp = names.index("bmi"), names.index("bp")
    X2 = X[:, [idx_bmi, idx_bp]]
    poly_names, sk_coefs, sk_intercept, ne_coefs, ne_intercept, coef_diff, intercept_diff = (
        r.polynomial_matches_normal_equations(X2, y, degree=2, feature_names=["bmi", "bp"])
    )
    print(f"  design matrix columns   : {poly_names}")
    print(f"  sklearn coefficients    : {sk_coefs}")
    print(f"  normal-eq coefficients  : {ne_coefs}")
    print(f"  sklearn intercept       : {sk_intercept}")
    print(f"  normal-eq intercept     : {ne_intercept}")
    print(f"  max abs coefficient gap : {coef_diff:.2e}")
    print(f"  max abs intercept gap   : {intercept_diff:.2e}")
    r2_with, r2_without, interaction_coef = r.interaction_term_effect(X2, y)
    print(f"  R2 with 'bmi bp' term   : {r2_with:.6f}")
    print(f"  R2 without it           : {r2_without:.6f}")
    print(f"  interaction coefficient : {interaction_coef:.6f}")

    rule("7. R2 never decreases when you add a predictor -- even noise")
    rows = r.r2_with_added_noise_columns(X, y, [1, 2, 5, 10], seed=42)
    print("  noise columns   R2         delta vs 10 real predictors")
    for n_noise, r2, delta in rows:
        print(f"  {n_noise:>13d}   {r2:.6f}   {delta:+.6f}")

    rule("8. Standardising changes the coefficients, not the model")
    raw_coefs, scaled_coefs, r2_raw, r2_scaled, max_pred_diff = r.scaling_effect(X, y)
    print("  name    raw coef     scaled coef")
    for name, raw_c, scaled_c in zip(names, raw_coefs, scaled_coefs):
        print(f"  {name:<5}   {raw_c:>9.4f}   {scaled_c:>10.4f}")
    print(f"  R2 raw / scaled                : {r2_raw:.6f} / {r2_scaled:.6f}")
    print(f"  max abs prediction difference  : {max_pred_diff:.2e}")


if __name__ == "__main__":
    main()

"""The reference solutions: what changes once a second predictor joins the first.

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

import numpy as np
import pytest

import regression_lib as r


@pytest.fixture(scope="module")
def diabetes():
    return r.load_raw_diabetes()


@pytest.fixture(scope="module")
def bmi_bp(diabetes):
    X, y, names = diabetes
    idx_bmi, idx_bp = names.index("bmi"), names.index("bp")
    return X[:, [idx_bmi, idx_bp]], y


# --- 1. Correlation and variance inflation --------------------------------


def test_01_variance_inflation_factors_flag_the_correlated_serum_measurements(diabetes):
    X, y, names = diabetes
    vifs = r.variance_inflation_factors(X, names)
    assert vifs == {
        "age": 1.2173,
        "sex": 1.2781,
        "bmi": 1.5094,
        "bp": 1.4594,
        "s1": 59.2025,
        "s2": 39.1934,
        "s3": 15.4022,
        "s4": 8.891,
        "s5": 10.076,
        "s6": 1.4846,
    }
    # Every clinical measurement (age, sex, bmi, bp) sits near 1: barely
    # explained by the other nine. Every serum measurement (s1-s5) sits
    # well above the common rule-of-thumb cutoff of 5.
    for name in ("age", "sex", "bmi", "bp"):
        assert vifs[name] < 2.0
    for name in ("s1", "s2", "s3", "s4", "s5"):
        assert vifs[name] > 5.0


def test_01b_s1_and_s2_are_the_most_correlated_predictor_pair(diabetes):
    X, y, names = diabetes
    assert round(r.correlation(X, names, "s1", "s2"), 4) == 0.8967
    assert round(r.correlation(X, names, "s3", "s4"), 4) == -0.7385
    # bmi and bp -- the two clinical measurements -- are nowhere near as
    # entangled with each other.
    assert abs(r.correlation(X, names, "bmi", "bp")) < 0.4


# --- 2. The centrepiece: an exact duplicate column -------------------------


def test_02_an_exact_duplicate_splits_the_coefficient_but_not_the_sum(diabetes):
    X, y, names = diabetes
    idx_s1 = names.index("s1")
    original, coef_a, coef_b, max_diff, r2_orig, r2_dup = r.duplicate_column_exact(X, y, idx_s1)
    assert round(original, 4) == -1.09
    assert round(coef_a, 4) == -0.545
    assert round(coef_b, 4) == -0.545
    # Neither half equals the original coefficient -- but their sum does,
    # to eight decimal places. The normal equations only ever "see" the
    # combined effect of two identical columns.
    assert abs((coef_a + coef_b) - original) < 1e-8


def test_02b_the_exact_duplicate_changes_nothing_about_the_model_itself(diabetes):
    X, y, names = diabetes
    idx_s1 = names.index("s1")
    _original, _a, _b, max_diff, r2_orig, r2_dup = r.duplicate_column_exact(X, y, idx_s1)
    # Every prediction the model makes is unchanged to eleven decimal places.
    assert max_diff < 1e-10
    assert abs(r2_dup - r2_orig) < 1e-10
    assert round(r2_orig, 4) == 0.5177


# --- 3. Breaking the tie with noise makes the split arbitrary --------------


def test_03_a_tiny_amount_of_noise_lets_the_two_coefficients_swing_wildly(diabetes):
    X, y, names = diabetes
    idx_s1 = names.index("s1")
    noise_scale = 0.01 * float(X[:, idx_s1].std())
    coef_a, coef_b, coef_sum, max_diff, r2 = r.duplicate_column_noisy(X, y, idx_s1, noise_scale, seed=0)
    assert round(coef_a, 4) == 0.7592
    assert round(coef_b, 4) == -1.8451
    # A one-percent noise perturbation was enough to send the original
    # -1.09 coefficient POSITIVE. The near-duplicate is no longer tied
    # exactly, so least squares picks a definite -- but essentially
    # arbitrary -- way to split the shared effect.
    assert coef_a > 0
    assert round(coef_sum, 4) == -1.0859
    assert round(r2, 4) == 0.5178


def test_03b_across_many_noise_draws_the_sum_and_the_predictions_hold_steady(diabetes):
    X, y, names = diabetes
    idx_s1 = names.index("s1")
    noise_scale = 0.01 * float(X[:, idx_s1].std())
    result = r.duplicate_noisy_spread(X, y, idx_s1, noise_scale, range(10))
    # The two individual coefficients range over more than twelve units --
    # further apart than the original coefficient is from zero, and they
    # cross zero repeatedly.
    assert result["coef_a"]["sd"] > 4.0
    assert result["coef_b"]["sd"] > 4.0
    assert result["coef_a"]["min"] < 0 < result["coef_a"]["max"]
    # Their sum barely moves: two orders of magnitude steadier than either
    # coefficient alone.
    assert result["sum"]["sd"] < 0.05
    assert round(result["sum"]["mean"], 2) == -1.09
    # And the predictions themselves move by a few units on a target whose
    # own spread is 77 -- noticeable, but nowhere near what the coefficient
    # swings would suggest.
    assert result["max_pred_diff_overall"] < 10.0
    assert result["r2"]["sd"] < 0.001


# --- 4. Instability tracks the VIF, not just the anecdote ------------------


def test_04_bootstrap_resampling_shows_high_vif_predictors_wobble_more(diabetes):
    X, y, names = diabetes
    boot = r.bootstrap_coefficient_spread(X, y, names, reps=500, seed=0)
    high_vif = ["s1", "s2", "s3", "s4"]
    low_vif = ["bmi", "bp", "sex"]
    high_cv = np.mean([boot[name]["cv"] for name in high_vif])
    low_cv = np.mean([boot[name]["cv"] for name in low_vif])
    # age's own coefficient of variation is enormous (4.70) because its
    # mean sits near zero, which inflates the ratio rather than reflecting
    # genuine instability -- excluded from this comparison for that reason.
    assert high_cv > low_cv
    assert boot["s1"]["cv"] > boot["bmi"]["cv"]
    assert boot["s3"]["cv"] > boot["bp"]["cv"]
    assert round(boot["s1"]["cv"], 2) == 0.51
    assert round(boot["bmi"]["cv"], 2) == 0.13


# --- 5. Holding the others constant can flip a sign -------------------------


def test_05_conditioning_on_the_other_nine_predictors_flips_four_signs(diabetes):
    X, y, names = diabetes
    result = r.simple_vs_multiple_coefficients(X, y, names)
    flips = {name for name, values in result.items() if values["sign_flip"]}
    assert flips == {"age", "sex", "s1", "s3"}
    # s1 is the headline: positive alone, negative once s2 is in the model.
    assert result["s1"]["simple"] == 0.4723
    assert result["s1"]["multiple"] == -1.09
    # bmi and s5 do not flip -- their relationship with the target survives
    # conditioning on the other nine predictors.
    assert result["bmi"]["sign_flip"] is False
    assert result["s5"]["sign_flip"] is False


# --- 6. A polynomial fit is linear in its parameters ------------------------


def test_06_polynomialfeatures_plus_linear_regression_matches_the_normal_equations(bmi_bp):
    X2, y = bmi_bp
    names, sk_coefs, sk_intercept, ne_coefs, ne_intercept, coef_diff, intercept_diff = (
        r.polynomial_matches_normal_equations(X2, y, degree=2, feature_names=["bmi", "bp"])
    )
    assert names == ["bmi", "bp", "bmi^2", "bmi bp", "bp^2"]
    assert sk_coefs == ne_coefs
    assert sk_intercept == ne_intercept
    # Two different solution methods for the identical linear system agree
    # to well beyond floating-point noise.
    assert coef_diff < 1e-9
    assert intercept_diff < 1e-9


def test_06b_dropping_the_interaction_term_costs_real_r_squared(bmi_bp):
    X2, y = bmi_bp
    r2_with, r2_without, interaction_coef = r.interaction_term_effect(X2, y)
    assert round(r2_with, 6) == 0.40417
    assert round(r2_without, 6) == 0.399896
    assert round(interaction_coef, 6) == 0.095079
    # bmi^2 and bp^2 describe each predictor curving on its own; only the
    # interaction term describes bmi's effect changing with bp's level --
    # dropping it is a real, measurable loss of fit, not a bookkeeping one.
    assert r2_with > r2_without


# --- 7. R-squared never decreases, even for pure noise ----------------------


def test_07_r_squared_never_decreases_when_you_add_a_predictor_even_noise(diabetes):
    X, y, names = diabetes
    rows = r.r2_with_added_noise_columns(X, y, [1, 2, 5, 10], seed=42)
    assert rows == [
        (1, 0.518064, 0.000316),
        (2, 0.523041, 0.005293),
        (5, 0.527615, 0.009867),
        (10, 0.532455, 0.014707),
    ]
    r2_values = [row[1] for row in rows]
    assert all(a < b for a, b in zip(r2_values, r2_values[1:]))
    # Ten columns of pure numpy noise, with no relationship to the target
    # whatsoever, bought 1.47 points of R2 for free.
    assert rows[-1][2] > 0.01


# --- 8. Scaling changes the coefficients, not the model ---------------------


def test_08_standardizing_changes_the_coefficients_not_the_predictions(diabetes):
    X, y, names = diabetes
    raw_coefs, scaled_coefs, r2_raw, r2_scaled, max_pred_diff = r.scaling_effect(X, y)
    # The two coefficient vectors are nowhere near each other in scale --
    # s1's raw coefficient is -1.09; its scaled coefficient is -37.68.
    idx_s1 = names.index("s1")
    assert raw_coefs[idx_s1] == -1.09
    assert scaled_coefs[idx_s1] == -37.68
    assert max(abs(c) for c in raw_coefs) < 70
    assert max(abs(c) for c in scaled_coefs) < 40
    # But the fit itself -- what it predicts, and how well -- is identical.
    assert r2_raw == r2_scaled
    assert max_pred_diff < 1e-9

"""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 regression_lib as r


def test_the_dataset_loads_with_ten_raw_unit_predictors():
    X, y, names = r.load_raw_diabetes()
    assert X.shape == (442, 10)
    assert y.shape == (442,)
    assert names == ["age", "sex", "bmi", "bp", "s1", "s2", "s3", "s4", "s5", "s6"]
    # scaled=False: real units, not the default mean-centred unit-norm columns.
    assert X[:, names.index("age")].min() >= 19
    assert X[:, names.index("age")].max() <= 79
    assert y.min() >= 25 and y.max() <= 346


def test_a_perfectly_uncorrelated_pair_has_a_vif_near_one():
    rng = np.random.default_rng(0)
    X = rng.normal(size=(500, 3))
    y = X[:, 0] * 2 + rng.normal(size=500) * 0.1
    vifs = r.variance_inflation_factors(X, ["a", "b", "c"])
    for name in ("a", "b", "c"):
        assert 0.9 < vifs[name] < 1.15


def test_an_exact_duplicate_of_a_column_has_infinite_vif():
    rng = np.random.default_rng(0)
    X = rng.normal(size=(200, 2))
    X_dup = np.hstack([X, X[:, [0]]])
    vifs = r.variance_inflation_factors(X_dup, ["a", "b", "a_copy"])
    assert vifs["a"] == float("inf")
    assert vifs["a_copy"] == float("inf")


def test_duplicating_an_unrelated_random_column_barely_moves_anything():
    rng = np.random.default_rng(1)
    X = rng.normal(size=(300, 2))
    y = X[:, 0] * 3.0 + rng.normal(size=300) * 0.5
    original_coef, coef_a, coef_b, max_diff, r2_orig, r2_dup = r.duplicate_column_exact(X, y, 1)
    # Column 1 has nothing to do with y, so splitting its (near-zero) effect
    # in half still leaves both halves small, and nothing else moves.
    assert abs(coef_a + coef_b - original_coef) < 1e-8
    assert max_diff < 1e-8
    assert abs(r2_dup - r2_orig) < 1e-8


def test_the_spread_helper_reports_mean_sd_min_max():
    result = r.spread([1.0, 2.0, 3.0, 4.0, 5.0])
    assert result == {"mean": 3.0, "sd": 1.4142, "min": 1.0, "max": 5.0}

lesson_id: D150
day: 150
kind: guided-build
languages:
  - python
  - bash
setup_commands:
  - cd labs/sections/machine-learning/day-150-multiple-and-polynomial-regression
  - python3 -m venv .venv
  - .venv/bin/pip install -r requirements/requirements.txt
  - >-
    .venv/bin/python3 -c "import numpy, sklearn; print(numpy.__version__,
    sklearn.__version__)"
run_commands:
  - .venv/bin/pytest examples -q
  - .venv/bin/pytest starter -q
  - .venv/bin/python3 examples/report_measurements.py
test_commands:
  - bash tests/run_tests.sh
cleanup_commands:
  - >-
    find . -path ./.venv -prune -o -type d -name '__pycache__' -print -exec rm -rf -- {}
    +
  - rm -rf .pytest_cache
  - 'rm -rf .venv  # optional: removes the lab virtual environment'
  - 'git checkout -- starter/  # optional: reset your work'
requires_network: true
requires_api_key: false
estimated_minutes: 60
last_executed: '2026-08-27'
executed_on: >-
  macOS 26.5.2 (Apple Silicon, arm64, CPU only -- no GPU is needed or used), Python
  3.14.0, numpy 2.5.2, scikit-learn 1.9.0, pytest 9.1.1, bash 3.2.57 -- bash
  tests/run_tests.sh -> 14 checks, 0 failure(s), exit 0. pytest examples -q -> 17 passed.
  pytest starter -q -> 5 passed, 12 skipped (the five machinery checks in
  test_regression_lib.py are solved in both directories; the twelve 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
  1.18.1, joblib 1.5.3 and threadpoolctl 3.6.0 as its own dependencies, none of which
  this lab imports directly. The lab is fully offline after the pip install -- the only
  dataset used, sklearn.datasets.load_diabetes(scaled=False), ships bundled inside the
  scikit-learn package as a compressed CSV and is never downloaded at run time, and
  harness check 9 confirms no URL appears anywhere in starter/ or examples/ source.
  Section 7 of the harness copies examples/ into a mktemp-d scratch directory, confirms
  17 passed, rewrites `assert round(interaction_coef, 6) == 0.095079` to 0.0, confirms a
  non-zero exit naming the failing test (test_06b_dropping_the_interaction_term_costs_real_r_squared),
  and removes the scratch directory. Separately, by hand,
  `assert vifs[name] > 5.0` for s1-s5 in test_01 was changed to `> 500.0` and the whole
  harness re-run: it reported 14 checks, 2 failure(s) and exited 1 (both the pytest run
  and the pytest-free direct reproduction in section 2 caught it, since s3's VIF of
  15.4022 is well under 500); the file was restored and the harness returned to 14
  checks, 0 failure(s), exit 0. MEASURED PAIRS, all captured verbatim in
  expected-output/measured-values.txt. (1) VARIANCE INFLATION, ten raw-unit predictors
  from load_diabetes(scaled=False): age 1.2173, sex 1.2781, bmi 1.5094, bp 1.4594, s1
  59.2025, s2 39.1934, s3 15.4022, s4 8.891, s5 10.076, s6 1.4846 -- every clinical
  measurement under 2, every one of s1-s5 over the rule-of-thumb cutoff of 5. Correlation
  s1-s2 is 0.8967, s3-s4 is -0.7385, bmi-bp is only 0.3954. (2) THE CENTREPIECE: appending
  an exact copy of s1 (original coefficient -1.09) splits it into two coefficients of
  -0.545 each -- their sum equals the original to eight decimal places, the maximum
  absolute prediction change is 3.98e-12, and R2 is unchanged at 0.5177. Breaking the
  exact tie with 1 percent Gaussian noise on the copy (noise scale 0.3457, on a target
  whose own standard deviation is 77) sends the split coefficients to +0.7592 and -1.8451
  at seed 0 -- the ORIGINAL was negative, one half is now positive -- while their sum
  stays at -1.0859 and R2 at 0.5178. Across ten noise seeds: both coefficients have
  standard deviation above 4.4 and cross zero (roughly -7 to +6), their sum has standard
  deviation 0.0144 and mean -1.0891, the largest single prediction move across all ten
  seeds is 6.5911, and R2's standard deviation is 0.0003. Wild coefficients, stable
  predictions. (3) BOOTSTRAP INSTABILITY, 500 resamples per predictor: coefficient of
  variation is 0.51 for s1 (VIF 59.2) against 0.13 for bmi (VIF 1.51) and 0.20 for bp
  (VIF 1.46) -- higher VIF predicts more bootstrap wobble. (4) SIGN FLIPS: conditioning
  on the other nine predictors flips the sign of age, sex, s1 and s3. s1 alone is +0.4723;
  s1 with the other nine held constant is -1.09. bmi, bp, s2, s4, s5, s6 do not flip.
  (5) POLYNOMIAL MECHANICS: PolynomialFeatures(degree=2) plus LinearRegression on bmi and
  bp gives coefficients [-1.891873, -2.116267, 0.024891, 0.095079, 0.004919] and intercept
  99.877329, identical to solving the normal equations on the same expanded design matrix
  by hand (max coefficient gap 1.49e-13, max intercept gap 9.81e-13). Dropping the
  interaction term (bmi bp) drops R2 from 0.404170 to 0.399896. (6) R2 NEVER DECREASES:
  adding 1, 2, 5 and 10 columns of pure numpy noise to the ten real predictors raises R2
  from a 10-predictor baseline of 0.5177 to 0.518064, 0.523041, 0.527615 and 0.532455 --
  strictly increasing every time, on columns with zero relationship to the target.
  (7) SCALING: standardising every predictor changes s1's coefficient from -1.09 to
  -37.68 -- more than 30 times larger in magnitude, and several other coefficients move
  by similar factors -- while R2 (0.517748) and every prediction (max difference
  1.14e-12) are unchanged. TWO HONESTY CALLS. FIRST: the noisy-duplicate result depends
  on which noise seed is quoted for its individual numbers (seed 0 sends one coefficient
  positive; other seeds send different individual values), so exercise 3 quotes one
  concrete seed to make the phenomenon vivid and exercise 3b reports the ten-seed spread
  -- mean, standard deviation, minimum and maximum -- as the defensible claim, following
  the same discipline Day 144 used for its temporal-split result. SECOND: age's bootstrap
  coefficient of variation (4.70) is the largest of any predictor in the dataset despite
  age having one of the lowest VIFs (1.22) -- because its own mean coefficient sits near
  zero, which inflates a sd-over-mean ratio without reflecting genuine instability. The
  VIF-instability comparison in exercise 4 excludes age for exactly this reason and says
  so in the assertion's own comment, rather than silently cherry-picking predictors that
  make the correlation look cleaner than it is. Harness check 8 re-runs the duplicate-
  column instability on a second predictor (s2) and at different noise seeds, the R2
  monotonicity at three different noise seeds, and the VIF-instability link at a
  different bootstrap replication count, so no directional claim in the lesson rests on
  the single seed it quotes.

# Requirements

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

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

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

## Why the versions are pinned exactly

`sklearn.linear_model.LinearRegression` solves the normal equations with a
deterministic least-squares routine, so most of this lab's numbers are
exact arithmetic rather than sampled draws -- but two things do depend on
the pins: `numpy.random.default_rng`, used for the noise columns and the
bootstrap resamples, carries no stream-compatibility guarantee across
NumPy versions, and scikit-learn's own internals (LAPACK routine choice,
convergence tolerances) can shift a coefficient in its last few decimal
places between minor versions.

What does not depend on the pins: every formula (variance inflation as
`1 / (1 - R2)`, the standard-error-free identity that two duplicate
coefficients sum to the original), every structural fact (the exact
duplicate's predictions matching to floating-point precision, the
polynomial fit matching the normal equations), and the direction of every
result -- a duplicated correlated predictor destabilises its coefficients
while leaving predictions almost untouched, high-VIF predictors wobble
more under resampling than low-VIF ones, R2 never decreases when a
predictor is added, and standardising changes coefficients without
changing predictions. Harness check 8 re-runs several of these directions
at seeds and predictors the lesson does not quote.

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

## Installing

From the lab directory:

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

The install step needs the network. Everything after it is offline: the
only dataset this lab uses is `sklearn.datasets.load_diabetes`, which
ships bundled inside the scikit-learn package you just installed --
nothing is downloaded at run time.

## 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 150 lab brief — Many Predictors, One Model

Day 148 gave you a line through one predictor. Day 149 gave you a reason
to square the error before minimising it. Neither told you what changes
once a second predictor joins the first -- and the honest answer is: more
than you would guess.

This lab measures it, on `sklearn.datasets.load_diabetes(scaled=False)`:
442 patients, ten predictors in their real units -- age in years, sex
coded 1 or 2, bmi, average blood pressure, and six serum measurements
`s1`-`s6` -- and one target, a quantitative measure of disease progression
a year after baseline.

## The claim you are here to measure

> A coefficient in a multiple regression means "the change in the target
> for one unit of this predictor, holding every other predictor fixed" --
> and that phrase is doing enormous work.

Two predictors that are correlated with each other, not just with the
target, can be held "fixed" relative to each other only approximately,
because moving one tends to move the other in real data. The regression
still fits -- the predictions can be excellent -- but the individual
coefficients become an accounting exercise between correlated columns
rather than a stable description of anything.

## The centrepiece, and the number to watch

`s1` and `s2`, two of the six serum measurements, correlate at **0.8967**.
Append an exact copy of `s1` to the design matrix and refit:

| | original | duplicate model |
| --- | --- | --- |
| `s1` coefficient | −1.0900 | −0.5450 |
| copy's coefficient | -- | −0.5450 |
| **sum** | −1.0900 | **−1.0900** |
| R2 | 0.5177 | 0.5177 |
| max prediction change | -- | 3.98 × 10⁻¹² |

Neither half matches the original coefficient. Their **sum** does, to
eight decimal places, because the normal equations only ever see the
*combined* effect of two identical columns and have no way to prefer one
split over another.

Now break the exact tie with one percent of noise and refit at ten
different noise seeds:

```text
  coefficient a : sd 4.4258, range -6.9523 to 5.6258
  coefficient b : sd 4.4141, range -6.6988 to 5.8570
  their sum     : sd 0.0144, mean -1.0891
  largest single prediction move across all ten seeds: 6.5911
```

**Wild coefficients, stable predictions.** That contrast is the whole
lesson. A model can be excellent at what it predicts and worthless as a
description of "the effect of `s1`" at the same time, and nothing about
its accuracy will tell you so.

## Variance inflation, computed directly

Exercise 1 computes the standard diagnostic: regress each predictor on
the other nine, and take `1 / (1 - R2)`.

| predictor | VIF | predictor | VIF |
| --- | --- | --- | --- |
| age | 1.2173 | s1 | 59.2025 |
| sex | 1.2781 | s2 | 39.1934 |
| bmi | 1.5094 | s3 | 15.4022 |
| bp | 1.4594 | s4 | 8.8910 |
| s6 | 1.4846 | s5 | 10.0760 |

The four clinical measurements sit near 1 -- barely explained by the other
nine. Every serum measurement sits above the common rule-of-thumb cutoff
of 5, and `s1` at 59.2 is the worst of them: 59 times the variance a truly
independent predictor would have.

## Two things that are true and easy to misread

**A polynomial fit is still a linear model.** "Linear regression" does not
mean "a straight line" -- it means linear *in the parameters*. Fit
`PolynomialFeatures(degree=2)` followed by `LinearRegression` on `bmi` and
`bp`, and solve the identical expanded design matrix by hand with
`numpy.linalg.lstsq`. The two answers agree to thirteen decimal places,
because they are solving the same linear system by two different routes.

**R2 never decreases when you add a predictor -- even a column of pure
noise.** Add 1, 2, 5 and 10 columns of `numpy.random.default_rng` noise
with no relationship to the target whatsoever, and R2 climbs from 0.5177
to 0.5325. Nothing was learned; the model simply has more knobs to turn.

## 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 twelve skips.
3. Replace one `pytest.skip(...)` at a time with real code. The skip text
   names the exact helper and the exact value to assert.
4. Print the measured pair in every exercise. A number you did not print
   is a number you did not look at.
5. When you want the whole measured table at once, run
   `.venv/bin/python3 examples/report_measurements.py`.

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

## What this lab deliberately does not cover

Ridge and lasso -- the standard fix for the instability you are about to
measure -- belong to Day 151. Adjusted R2, the fix for the never-decreases
problem in exercise 7, belongs to Day 152. Neither is implemented here;
this lab only measures the two problems those days solve.

"""Many predictors, measured: what changes when there is more than one.

Day 148 covered the line through one predictor. Day 149 covered why the
loss is squared error. This module measures what is new once a second
predictor joins the first: a coefficient's meaning becomes conditional on
"holding the others constant", and when predictors are correlated with
each other -- not just with the target -- that condition gets expensive.

The centrepiece is a duplicated predictor. Two columns that carry the same
information cannot be told apart by the normal equations, so the fit
smears one true effect across both coefficients in whatever proportion
happens to minimise squared error that day -- while the *sum* of the two
coefficients, and every prediction the model makes, barely moves at all.

Everything here is deterministic given a seed, and every dataset is either
the bundled `sklearn.datasets.load_diabetes(scaled=False)` -- ten raw-unit
clinical predictors, real measurements, no download -- or built from it.
"""

from __future__ import annotations

import numpy as np

from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures, StandardScaler


# --------------------------------------------------------------------------
# 0. The dataset itself
# --------------------------------------------------------------------------


def load_raw_diabetes():
    """The ten raw-unit predictors and the target, with feature names.

    `scaled=False` matters: it returns age in years, sex coded 1/2, bmi,
    average blood pressure, and six serum measurements s1-s6 in their
    original units, which is what makes a coefficient's magnitude mean
    anything. The default `scaled=True` mean-centres and unit-norm scales
    every column before you ever see it -- exercise 8 measures exactly
    what that substitution does and does not change.
    """
    bunch = load_diabetes(scaled=False)
    return bunch.data, bunch.target, list(bunch.feature_names)


def fit(X, y):
    """A plain ordinary-least-squares fit -- the one tool this lesson uses."""
    return LinearRegression().fit(X, y)


# --------------------------------------------------------------------------
# 1. Correlation and variance inflation
# --------------------------------------------------------------------------


def correlation(X, names, a: str, b: str) -> float:
    """The Pearson correlation between two named predictor columns."""
    i, j = names.index(a), names.index(b)
    return float(np.corrcoef(X[:, i], X[:, j])[0, 1])


def variance_inflation_factors(X, names) -> dict:
    """VIF for every column: regress it on the rest, VIF = 1 / (1 - R2).

    A VIF of 1 means a predictor is unrelated to the others. Above 5 or 10
    is the usual rule of thumb for "correlated enough to worry about" --
    this function does not apply the rule, it only computes the number the
    rule is applied to.
    """
    result = {}
    for i, name in enumerate(names):
        others = np.delete(X, i, axis=1)
        target = X[:, i]
        r2 = LinearRegression().fit(others, target).score(others, target)
        result[name] = float("inf") if r2 >= 1.0 else round(1.0 / (1.0 - r2), 4)
    return result


# --------------------------------------------------------------------------
# 2. The centrepiece: an exact duplicate column
# --------------------------------------------------------------------------


def duplicate_column_exact(X, y, col_index: int):
    """Append an exact copy of one column and refit.

    Returns ``(original_coef, dup_coef_a, dup_coef_b, max_abs_pred_diff,
    r2_original, r2_dup)``. The two duplicate coefficients need not equal
    the original one individually -- only their sum does, because the
    normal equations only ever "see" the combined effect of two identical
    columns, and there is no unique way to split it.
    """
    original = fit(X, y)
    X_dup = np.hstack([X, X[:, [col_index]]])
    dup = fit(X_dup, y)
    pred_original = original.predict(X)
    pred_dup = dup.predict(X_dup)
    return (
        float(original.coef_[col_index]),
        float(dup.coef_[col_index]),
        float(dup.coef_[-1]),
        float(np.max(np.abs(pred_dup - pred_original))),
        float(original.score(X, y)),
        float(dup.score(X_dup, y)),
    )


def duplicate_column_noisy(X, y, col_index: int, noise_scale: float, seed: int):
    """Append a near-duplicate: the same column plus a little noise.

    Breaking the exact tie lets ordinary least squares pick a *unique*
    split of the combined effect again -- but which split it picks depends
    on which way the noise happened to fall, which is the whole point.
    Returns ``(coef_a, coef_b, sum_coefs, max_abs_pred_diff, r2)``.
    """
    original = fit(X, y)
    pred_original = original.predict(X)
    rng = np.random.default_rng(seed)
    noise = rng.normal(scale=noise_scale, size=X.shape[0])
    near_dup = X[:, col_index] + noise
    X_noisy = np.hstack([X, near_dup.reshape(-1, 1)])
    model = fit(X_noisy, y)
    pred = model.predict(X_noisy)
    return (
        float(model.coef_[col_index]),
        float(model.coef_[-1]),
        float(model.coef_[col_index] + model.coef_[-1]),
        float(np.max(np.abs(pred - pred_original))),
        float(model.score(X_noisy, y)),
    )


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


def duplicate_noisy_spread(X, y, col_index: int, noise_scale: float, seeds) -> dict:
    """Across many noise draws: how the two coefficients wander, and how
    their sum, the predictions and R2 do not.

    Returns a dict with ``coef_a``, ``coef_b`` and ``sum`` spreads (each
    from :func:`spread`), plus ``max_pred_diff`` (the largest single
    prediction movement seen across every seed) and an ``r2`` spread.
    """
    coef_a, coef_b, coef_sum, max_diffs, r2s = [], [], [], [], []
    for seed in seeds:
        a, b, total, max_diff, r2 = duplicate_column_noisy(X, y, col_index, noise_scale, seed)
        coef_a.append(a)
        coef_b.append(b)
        coef_sum.append(total)
        max_diffs.append(max_diff)
        r2s.append(r2)
    return {
        "coef_a": spread(coef_a),
        "coef_b": spread(coef_b),
        "sum": spread(coef_sum),
        "max_pred_diff_overall": round(max(max_diffs), 4),
        "r2": spread(r2s),
    }


# --------------------------------------------------------------------------
# 3. Bootstrap coefficient instability, across all ten predictors
# --------------------------------------------------------------------------


def bootstrap_coefficient_spread(X, y, names, reps: int = 500, seed: int = 0) -> dict:
    """Resample the rows with replacement, refit, and record each coefficient.

    Returns ``{name: {"mean": ..., "sd": ..., "cv": ...}}`` where ``cv`` is
    the coefficient of variation, ``abs(sd / mean)`` -- a scale-free way to
    compare how much a coefficient wanders relative to its own size.
    """
    rng = np.random.default_rng(seed)
    n = X.shape[0]
    collected = {name: [] for name in names}
    for _ in range(reps):
        idx = rng.integers(0, n, size=n)
        model = fit(X[idx], y[idx])
        for i, name in enumerate(names):
            collected[name].append(model.coef_[i])
    result = {}
    for name in names:
        arr = np.asarray(collected[name])
        mean = float(arr.mean())
        sd = float(arr.std())
        result[name] = {
            "mean": round(mean, 4),
            "sd": round(sd, 4),
            "cv": round(abs(sd / mean), 4) if mean != 0 else None,
        }
    return result


# --------------------------------------------------------------------------
# 4. What "holding the others constant" changes: simple vs. multiple
# --------------------------------------------------------------------------


def simple_vs_multiple_coefficients(X, y, names) -> dict:
    """Each predictor's coefficient alone, and again inside the full model.

    Returns ``{name: {"simple": ..., "multiple": ..., "sign_flip": bool}}``.
    A sign flip means the predictor's *apparent* relationship with the
    target, taken alone, points the opposite way from its *conditional*
    relationship once the other nine predictors are held constant.
    """
    result = {}
    full = fit(X, y)
    for i, name in enumerate(names):
        alone = fit(X[:, [i]], y)
        simple_coef = float(alone.coef_[0])
        multiple_coef = float(full.coef_[i])
        result[name] = {
            "simple": round(simple_coef, 4),
            "multiple": round(multiple_coef, 4),
            "sign_flip": bool(np.sign(simple_coef) != np.sign(multiple_coef)),
        }
    return result


# --------------------------------------------------------------------------
# 5. A polynomial fit is linear in its parameters
# --------------------------------------------------------------------------


def polynomial_matches_normal_equations(X2, y, degree: int = 2, feature_names=None):
    """PolynomialFeatures + LinearRegression against a direct normal-equations solve.

    Builds the degree-``degree`` design matrix two ways: once through
    scikit-learn's transformer-plus-estimator pipeline, once by hand with
    ``numpy.linalg.lstsq`` on the identical expanded matrix. If a
    polynomial fit really is "linear in its parameters, not in x", the two
    must agree to floating-point precision, because they are solving the
    same linear system.

    Returns ``(feature_names, sklearn_coefs, sklearn_intercept,
    normal_eq_coefs, normal_eq_intercept, max_abs_coef_diff,
    max_abs_intercept_diff)``.
    """
    poly = PolynomialFeatures(degree=degree, include_bias=False)
    X_poly = poly.fit_transform(X2)
    expanded_names = list(poly.get_feature_names_out(feature_names))

    sklearn_model = fit(X_poly, y)

    design = np.hstack([np.ones((X_poly.shape[0], 1)), X_poly])
    beta, *_ = np.linalg.lstsq(design, y, rcond=None)

    coef_diff = float(np.max(np.abs(beta[1:] - sklearn_model.coef_)))
    intercept_diff = float(abs(beta[0] - sklearn_model.intercept_))

    return (
        expanded_names,
        [round(float(c), 6) for c in sklearn_model.coef_],
        round(float(sklearn_model.intercept_), 6),
        [round(float(b), 6) for b in beta[1:]],
        round(float(beta[0]), 6),
        coef_diff,
        intercept_diff,
    )


def interaction_term_effect(X2, y):
    """Compare a degree-2 fit with and without its interaction term.

    ``bmi^2`` and ``bp^2`` describe how each predictor curves on its own.
    ``bmi bp`` describes something neither can: whether the *effect* of
    one depends on the level of the other. Drop it and refit on the
    remaining four columns; the gap in R2 is what the interaction term was
    buying.

    Returns ``(r2_with_interaction, r2_without_interaction,
    interaction_coefficient)``.
    """
    poly = PolynomialFeatures(degree=2, include_bias=False)
    X_poly = poly.fit_transform(X2)
    feature_names = list(poly.get_feature_names_out(["a", "b"]))
    interaction_index = feature_names.index("a b")

    with_interaction = fit(X_poly, y)
    X_without = np.delete(X_poly, interaction_index, axis=1)
    without_interaction = fit(X_without, y)

    return (
        float(with_interaction.score(X_poly, y)),
        float(without_interaction.score(X_without, y)),
        float(with_interaction.coef_[interaction_index]),
    )


# --------------------------------------------------------------------------
# 6. R-squared never decreases when you add a predictor -- even noise
# --------------------------------------------------------------------------


def r2_with_added_noise_columns(X, y, noise_counts, seed: int = 42) -> list:
    """R2 of the fit as pure-noise columns are appended, one count at a time.

    Returns rows of ``(n_noise_columns, r2, delta_from_baseline)``. Every
    added column is `numpy.random.default_rng` noise with no relationship
    to the target whatsoever -- Day 152 owns the fix (adjusted R2); this
    function only measures the problem it fixes.
    """
    base_r2 = fit(X, y).score(X, y)
    rng = np.random.default_rng(seed)
    rows = []
    for n_noise in noise_counts:
        noise_columns = rng.normal(size=(X.shape[0], n_noise))
        X_aug = np.hstack([X, noise_columns])
        r2 = fit(X_aug, y).score(X_aug, y)
        rows.append((n_noise, round(r2, 6), round(r2 - base_r2, 6)))
    return rows


# --------------------------------------------------------------------------
# 7. Scaling changes the coefficients, not the model
# --------------------------------------------------------------------------


def scaling_effect(X, y):
    """Fit on raw units and on standardised units; compare everything.

    Returns ``(raw_coefs, scaled_coefs, raw_r2, scaled_r2,
    max_abs_pred_diff)``. Standardising centres and unit-variance-scales
    every column before fitting, which changes what one unit of a
    predictor means and therefore changes every coefficient's size -- but
    changes nothing about what the model actually predicts.
    """
    raw_model = fit(X, y)
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)
    scaled_model = fit(X_scaled, y)

    pred_raw = raw_model.predict(X)
    pred_scaled = scaled_model.predict(X_scaled)

    return (
        [round(float(c), 4) for c in raw_model.coef_],
        [round(float(c), 4) for c in scaled_model.coef_],
        float(raw_model.score(X, y)),
        float(scaled_model.score(X_scaled, y)),
        float(np.max(np.abs(pred_raw - pred_scaled))),
    )

"""Twelve exercises in what changes once a second predictor joins the first.

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 numpy as np  # noqa: F401  (you will need it)
import pytest

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


@pytest.fixture(scope="module")
def diabetes():
    return r.load_raw_diabetes()


@pytest.fixture(scope="module")
def bmi_bp(diabetes):
    X, y, names = diabetes
    idx_bmi, idx_bp = names.index("bmi"), names.index("bp")
    return X[:, [idx_bmi, idx_bp]], y


def test_01_variance_inflation_factors_flag_the_correlated_serum_measurements(diabetes):
    pytest.skip(
        "Assert r.variance_inflation_factors(X, names) equals the ten-entry "
        "dict in expected-output/measured-values.txt: age 1.2173, sex "
        "1.2781, bmi 1.5094, bp 1.4594, s1 59.2025, s2 39.1934, s3 15.4022, "
        "s4 8.891, s5 10.076, s6 1.4846. Then assert every one of age, sex, "
        "bmi, bp is under 2.0, and every one of s1-s5 is over 5.0 -- the "
        "usual VIF rule-of-thumb cutoff."
    )


def test_01b_s1_and_s2_are_the_most_correlated_predictor_pair(diabetes):
    pytest.skip(
        "Assert r.correlation(X, names, 's1', 's2') rounds to 0.8967 and "
        "r.correlation(X, names, 's3', 's4') rounds to -0.7385. Then assert "
        "abs(r.correlation(X, names, 'bmi', 'bp')) is under 0.4 -- the two "
        "clinical measurements are nowhere near as entangled as the serum "
        "measurements are with each other."
    )


def test_02_an_exact_duplicate_splits_the_coefficient_but_not_the_sum(diabetes):
    pytest.skip(
        "Call r.duplicate_column_exact(X, y, names.index('s1')). Assert the "
        "original coefficient rounds to -1.09 and the two duplicate "
        "coefficients each round to -0.545. Then assert their SUM equals "
        "the original coefficient to within 1e-8 -- neither half matches "
        "the original alone, but the combined effect is conserved exactly."
    )


def test_02b_the_exact_duplicate_changes_nothing_about_the_model_itself(diabetes):
    pytest.skip(
        "From the same call, assert the maximum absolute difference in "
        "predictions is under 1e-10 and the R2 values agree to within "
        "1e-10, with R2 rounding to 0.5177. Adding a column that carries no "
        "new information changes the model's arithmetic without changing "
        "what it predicts."
    )


def test_03_a_tiny_amount_of_noise_lets_the_two_coefficients_swing_wildly(diabetes):
    pytest.skip(
        "With noise_scale = 0.01 * X[:, s1_index].std(), call "
        "r.duplicate_column_noisy(X, y, s1_index, noise_scale, seed=0). "
        "Assert the first coefficient rounds to 0.7592 (POSITIVE -- the "
        "original was -1.09) and the second to -1.8451. Assert their sum "
        "rounds to -1.0859 and R2 rounds to 0.5178. Breaking the exact tie "
        "with a one-percent perturbation was enough to send the sign the "
        "wrong way."
    )


def test_03b_across_many_noise_draws_the_sum_and_the_predictions_hold_steady(diabetes):
    pytest.skip(
        "Call r.duplicate_noisy_spread(X, y, s1_index, noise_scale, "
        "range(10)). Assert both individual coefficients have a standard "
        "deviation over 4.0 and cross zero (min negative, max positive). "
        "Then assert the SUM's standard deviation is under 0.05 and its "
        "mean rounds to -1.09 at 2 decimals. Assert max_pred_diff_overall "
        "is under 10.0 and the R2 standard deviation is under 0.001. Wild "
        "coefficients, stable predictions -- that contrast is the lesson."
    )


def test_04_bootstrap_resampling_shows_high_vif_predictors_wobble_more(diabetes):
    pytest.skip(
        "Call r.bootstrap_coefficient_spread(X, y, names, reps=500, "
        "seed=0). Average the 'cv' (coefficient of variation) for "
        "['s1','s2','s3','s4'] and separately for ['bmi','bp','sex'], and "
        "assert the high-VIF average exceeds the low-VIF average. Then "
        "assert boot['s1']['cv'] rounds to 0.51 and boot['bmi']['cv'] "
        "rounds to 0.13 -- a predictor's VIF predicts how much its own "
        "coefficient wobbles under resampling."
    )


def test_05_conditioning_on_the_other_nine_predictors_flips_four_signs(diabetes):
    pytest.skip(
        "Call r.simple_vs_multiple_coefficients(X, y, names). Collect the "
        "names whose 'sign_flip' is True and assert the set equals "
        "{'age', 'sex', 's1', 's3'}. Then assert result['s1']['simple'] == "
        "0.4723 and result['s1']['multiple'] == -1.09 -- positive alone, "
        "negative once s2 (its correlated partner) is held constant."
    )


def test_06_polynomialfeatures_plus_linear_regression_matches_the_normal_equations(bmi_bp):
    pytest.skip(
        "Call r.polynomial_matches_normal_equations(X2, y, degree=2, "
        "feature_names=['bmi', 'bp']). Assert the feature names are "
        "['bmi', 'bp', 'bmi^2', 'bmi bp', 'bp^2'], the sklearn coefficients "
        "equal the normal-equations coefficients exactly (they were both "
        "rounded to 6 places), and both difference values are under 1e-9. "
        "A polynomial fit is linear in ITS parameters, and this is the "
        "proof: two different solvers on the identical design matrix agree."
    )


def test_06b_dropping_the_interaction_term_costs_real_r_squared(bmi_bp):
    pytest.skip(
        "Call r.interaction_term_effect(X2, y). Assert r2_with rounds to "
        "0.40417 and r2_without rounds to 0.399896, with r2_with strictly "
        "greater. Assert the interaction coefficient rounds to 0.095079. "
        "bmi^2 and bp^2 describe each predictor curving on its own; only "
        "'bmi bp' describes bmi's effect changing with bp's level."
    )


def test_07_r_squared_never_decreases_when_you_add_a_predictor_even_noise(diabetes):
    pytest.skip(
        "Call r.r2_with_added_noise_columns(X, y, [1, 2, 5, 10], seed=42) "
        "and assert it equals the four rows in "
        "expected-output/measured-values.txt, from (1, 0.518064, 0.000316) "
        "to (10, 0.532455, 0.014707). Assert the R2 column is strictly "
        "increasing. Every added column is pure numpy noise with no "
        "relationship to the target at all."
    )


def test_08_standardizing_changes_the_coefficients_not_the_predictions(diabetes):
    pytest.skip(
        "Call r.scaling_effect(X, y). Assert raw_coefs[s1_index] == -1.09 "
        "and scaled_coefs[s1_index] == -37.68 -- more than 30 times larger "
        "in magnitude. Then assert r2_raw == r2_scaled exactly and "
        "max_pred_diff is under 1e-9. Standardising changes what one unit "
        "of a predictor means, and therefore every coefficient's size -- "
        "and changes nothing about what the model predicts."
    )

"""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 regression_lib as r


def test_the_dataset_loads_with_ten_raw_unit_predictors():
    X, y, names = r.load_raw_diabetes()
    assert X.shape == (442, 10)
    assert y.shape == (442,)
    assert names == ["age", "sex", "bmi", "bp", "s1", "s2", "s3", "s4", "s5", "s6"]
    # scaled=False: real units, not the default mean-centred unit-norm columns.
    assert X[:, names.index("age")].min() >= 19
    assert X[:, names.index("age")].max() <= 79
    assert y.min() >= 25 and y.max() <= 346


def test_a_perfectly_uncorrelated_pair_has_a_vif_near_one():
    rng = np.random.default_rng(0)
    X = rng.normal(size=(500, 3))
    y = X[:, 0] * 2 + rng.normal(size=500) * 0.1
    vifs = r.variance_inflation_factors(X, ["a", "b", "c"])
    for name in ("a", "b", "c"):
        assert 0.9 < vifs[name] < 1.15


def test_an_exact_duplicate_of_a_column_has_infinite_vif():
    rng = np.random.default_rng(0)
    X = rng.normal(size=(200, 2))
    X_dup = np.hstack([X, X[:, [0]]])
    vifs = r.variance_inflation_factors(X_dup, ["a", "b", "a_copy"])
    assert vifs["a"] == float("inf")
    assert vifs["a_copy"] == float("inf")


def test_duplicating_an_unrelated_random_column_barely_moves_anything():
    rng = np.random.default_rng(1)
    X = rng.normal(size=(300, 2))
    y = X[:, 0] * 3.0 + rng.normal(size=300) * 0.5
    original_coef, coef_a, coef_b, max_diff, r2_orig, r2_dup = r.duplicate_column_exact(X, y, 1)
    # Column 1 has nothing to do with y, so splitting its (near-zero) effect
    # in half still leaves both halves small, and nothing else moves.
    assert abs(coef_a + coef_b - original_coef) < 1e-8
    assert max_diff < 1e-8
    assert abs(r2_dup - r2_orig) < 1e-8


def test_the_spread_helper_reports_mean_sd_min_max():
    result = r.spread([1.0, 2.0, 3.0, 4.0, 5.0])
    assert result == {"mean": 3.0, "sd": 1.4142, "min": 1.0, "max": 5.0}

#!/usr/bin/env bash
# Day 150 lab harness: "Many Predictors, One Model"
#
# 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

import regression_lib as r

errors = []


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


X, y, names = r.load_raw_diabetes()

# 1. VIF and correlation
vifs = r.variance_inflation_factors(X, names)
expect(
    "variance inflation factors",
    vifs,
    {
        "age": 1.2173, "sex": 1.2781, "bmi": 1.5094, "bp": 1.4594,
        "s1": 59.2025, "s2": 39.1934, "s3": 15.4022, "s4": 8.891,
        "s5": 10.076, "s6": 1.4846,
    },
)
expect("correlation s1,s2", round(r.correlation(X, names, "s1", "s2"), 4), 0.8967)
expect("correlation s3,s4", round(r.correlation(X, names, "s3", "s4"), 4), -0.7385)

# 2. Exact duplicate
idx_s1 = names.index("s1")
original, coef_a, coef_b, max_diff, r2_orig, r2_dup = r.duplicate_column_exact(X, y, idx_s1)
expect("original s1 coefficient", round(original, 4), -1.09)
expect("duplicate coefficient a", round(coef_a, 4), -0.545)
expect("duplicate coefficient b", round(coef_b, 4), -0.545)
if abs((coef_a + coef_b) - original) >= 1e-8:
    errors.append("the two duplicate coefficients did not sum to the original")
if max_diff >= 1e-10:
    errors.append(f"exact duplicate moved predictions by {max_diff}")
if abs(r2_dup - r2_orig) >= 1e-10:
    errors.append("exact duplicate changed R2")

# 3. Noisy duplicate
noise_scale = 0.01 * float(X[:, idx_s1].std())
n_coef_a, n_coef_b, n_sum, n_max_diff, n_r2 = r.duplicate_column_noisy(X, y, idx_s1, noise_scale, seed=0)
expect("noisy coefficient a (seed 0)", round(n_coef_a, 4), 0.7592)
expect("noisy coefficient b (seed 0)", round(n_coef_b, 4), -1.8451)
expect("noisy sum (seed 0)", round(n_sum, 4), -1.0859)
spread10 = r.duplicate_noisy_spread(X, y, idx_s1, noise_scale, range(10))
if spread10["coef_a"]["sd"] <= 4.0 or spread10["coef_b"]["sd"] <= 4.0:
    errors.append("noisy coefficients were less volatile than expected")
if spread10["sum"]["sd"] >= 0.05:
    errors.append(f"the sum was less stable than expected: sd={spread10['sum']['sd']}")
expect("sum mean, seeds 0-9 (2dp)", round(spread10["sum"]["mean"], 2), -1.09)
if spread10["max_pred_diff_overall"] >= 10.0:
    errors.append("predictions moved more than expected across noise seeds")

# 4. Bootstrap instability
boot = r.bootstrap_coefficient_spread(X, y, names, reps=500, seed=0)
expect("bootstrap s1 cv (2dp)", round(boot["s1"]["cv"], 2), 0.51)
expect("bootstrap bmi cv (2dp)", round(boot["bmi"]["cv"], 2), 0.13)
if boot["s1"]["cv"] <= boot["bmi"]["cv"]:
    errors.append("high-VIF predictor was not more unstable than a low-VIF one")

# 5. Sign flips
svm = r.simple_vs_multiple_coefficients(X, y, names)
flips = {name for name, v in svm.items() if v["sign_flip"]}
expect("predictors whose sign flips", flips, {"age", "sex", "s1", "s3"})
expect("s1 simple coefficient", svm["s1"]["simple"], 0.4723)
expect("s1 multiple coefficient", svm["s1"]["multiple"], -1.09)

# 6. Polynomial equals normal equations
idx_bmi, idx_bp = names.index("bmi"), names.index("bp")
X2 = X[:, [idx_bmi, idx_bp]]
poly_names, sk_coefs, sk_intercept, ne_coefs, ne_intercept, coef_diff, intercept_diff = (
    r.polynomial_matches_normal_equations(X2, y, degree=2, feature_names=["bmi", "bp"])
)
expect("expanded design matrix columns", poly_names, ["bmi", "bp", "bmi^2", "bmi bp", "bp^2"])
expect("sklearn coefficients == normal-eq coefficients", sk_coefs, ne_coefs)
if coef_diff >= 1e-9 or intercept_diff >= 1e-9:
    errors.append("sklearn and the normal equations disagreed beyond floating-point noise")

r2_with, r2_without, interaction_coef = r.interaction_term_effect(X2, y)
expect("R2 with interaction term", round(r2_with, 6), 0.40417)
expect("R2 without interaction term", round(r2_without, 6), 0.399896)
if r2_with <= r2_without:
    errors.append("the interaction term did not improve R2")

# 7. R2 with added noise columns
rows = r.r2_with_added_noise_columns(X, y, [1, 2, 5, 10], seed=42)
expect(
    "R2 with added noise columns",
    rows,
    [
        (1, 0.518064, 0.000316),
        (2, 0.523041, 0.005293),
        (5, 0.527615, 0.009867),
        (10, 0.532455, 0.014707),
    ],
)
r2_col = [row[1] for row in rows]
if not all(a < b for a, b in zip(r2_col, r2_col[1:])):
    errors.append("R2 did not strictly increase as noise columns were added")

# 8. Scaling
raw_coefs, scaled_coefs, r2_raw, r2_scaled, max_pred_diff = r.scaling_effect(X, y)
expect("raw s1 coefficient", raw_coefs[idx_s1], -1.09)
expect("scaled s1 coefficient", scaled_coefs[idx_s1], -37.68)
if r2_raw != r2_scaled:
    errors.append("scaling changed R2, which it must not")
if max_pred_diff >= 1e-9:
    errors.append("scaling changed predictions, which it must not")

if errors:
    for e in errors:
        print("ERROR:", e)
    sys.exit(1)
print("all direct checks passed")
PYEOF
)
if echo "$DIRECT_CHECK" | grep -q "all direct checks passed"; then
  ok "exercises 1-8 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 "^17 passed"; then
  ok "pytest examples -q -> 17 passed"
else
  fail "pytest examples -q did not report 17 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, 12 skipped"; then
  ok "pytest starter -q -> 5 passed, 12 skipped (the machinery checks pass; the twelve exercises are stubs)"
else
  fail "pytest starter -q did not report 5 passed, 12 skipped"
  echo "$STARTER_OUT" | tail -20 | sed 's/^/    /'
fi

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

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

echo ""
echo "7. Proof the harness can fail"
SCRATCH=$(mktemp -d "${TMPDIR:-/tmp}/d150-scratch.XXXXXX")
cp examples/*.py "$SCRATCH"/
SCRATCH_OUT=$("$PYTEST" "$SCRATCH" -q 2>&1)
if echo "$SCRATCH_OUT" | tail -1 | grep -qE "^17 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 round(interaction_coef, 6) == 0.095079'
replacement = 'assert round(interaction_coef, 6) == 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_06b_dropping_the_interaction_term_costs_real_r_squared"; then
  ok "breaking exercise 6b's assertion produces a non-zero exit and names the failing test"
else
  fail "broken copy did not fail as expected (exit=$BROKEN_STATUS)"
fi
rm -rf "$SCRATCH"

echo ""
echo "8. Key results hold at seeds and predictors the lesson does not quote"
DIRECTION=$("$PYTHON" - <<'PYEOF'
import sys
sys.path.insert(0, "examples")
import numpy as np
import regression_lib as r

problems = []
X, y, names = r.load_raw_diabetes()

# The duplicate-column contrast is not a property of one noise seed.
idx_s1 = names.index("s1")
noise_scale = 0.01 * float(X[:, idx_s1].std())
for seeds in ([10, 11, 12, 13, 14], [20, 21, 22, 23, 24]):
    spread = r.duplicate_noisy_spread(X, y, idx_s1, noise_scale, seeds)
    if spread["coef_a"]["sd"] <= 2.0:
        problems.append(f"seeds {seeds}: coefficient a was not volatile ({spread['coef_a']['sd']})")
    if spread["sum"]["sd"] >= 0.1:
        problems.append(f"seeds {seeds}: the sum was not stable ({spread['sum']['sd']})")

# The duplicate contrast is not specific to s1: try s2 (also high VIF).
idx_s2 = names.index("s2")
noise_scale_s2 = 0.01 * float(X[:, idx_s2].std())
spread_s2 = r.duplicate_noisy_spread(X, y, idx_s2, noise_scale_s2, range(5))
if spread_s2["coef_a"]["sd"] <= 1.0:
    problems.append("duplicating s2 did not produce unstable coefficients")
if spread_s2["sum"]["sd"] >= 0.1:
    problems.append("duplicating s2 did not leave the sum stable")

# R2 never decreasing is not a property of seed 42 alone.
for seed in (1, 2, 3):
    rows = r.r2_with_added_noise_columns(X, y, [1, 3, 8], seed=seed)
    values = [row[1] for row in rows]
    if not all(a < b for a, b in zip(values, values[1:])):
        problems.append(f"seed {seed}: R2 did not strictly increase with added noise columns")

# Bootstrap instability tracking VIF is not a property of one replication count.
boot = r.bootstrap_coefficient_spread(X, y, names, reps=150, seed=3)
if boot["s2"]["cv"] <= boot["bp"]["cv"]:
    problems.append("at a different replication count, a high-VIF predictor was not more unstable")

if problems:
    for p in problems:
        print("ERROR:", p)
else:
    print("every direction held")
PYEOF
)
if [ "$DIRECTION" = "every direction held" ]; then
  ok "duplicate-column instability, R2 monotonicity and VIF-linked instability hold beyond the quoted seeds and predictors"
else
  fail "a direction failed beyond the quoted seed"
  echo "$DIRECTION" | sed 's/^/    /'
fi

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

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