Hands-on lab — Day 151: Regularization: Ridge and Lasso

"""Ridge and lasso, measured: what the penalty does, and why the two differ
in kind rather than degree.

Day 145 already measured that a ridge penalty rescued an overfit polynomial
by a factor of 39,588, and that a penalty trades variance for bias. This
module does not re-measure that trade; it measures the CONTRAST between the
two penalty shapes: an L2 penalty that shrinks every coefficient toward zero
but never quite reaches it, and an L1 penalty that drives some coefficients
to exactly zero and leaves the rest alone.

Everything here is deterministic given a seed, and every dataset is either
scikit-learn's bundled ``load_diabetes`` or generated on the spot with
``numpy.random.default_rng`` or ``sklearn.datasets.make_regression``.
"""

from __future__ import annotations

import numpy as np

from sklearn.datasets import load_diabetes, make_regression
from sklearn.linear_model import ElasticNet, Lasso, LassoCV, Ridge
from sklearn.metrics import r2_score
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

FEATURE_NAMES = ["age", "sex", "bmi", "bp", "s1", "s2", "s3", "s4", "s5", "s6"]


# --------------------------------------------------------------------------
# 1. Ridge zeros nothing; lasso zeros progressively more
# --------------------------------------------------------------------------


def load_train_test(scaled: bool = True, seed: int = 0):
    """The diabetes dataset, split the way every exercise in this lab uses.

    ``scaled=True`` is scikit-learn's own bundled version: mean-centred and
    scaled so every column has unit L2 norm. ``scaled=False`` returns the raw
    measurement units (age in years, bmi as a ratio, bp in mm Hg, and so on).
    """
    X, y = load_diabetes(return_X_y=True, scaled=scaled)
    return train_test_split(X, y, test_size=0.25, random_state=seed)


def zero_counts_and_r2(alphas):
    """For each alpha, how many lasso coefficients are exactly zero, how
    many ridge coefficients are exactly zero, and each model's test R2.

    Returns rows of ``(alpha, lasso_zeros, lasso_r2, ridge_zeros, ridge_r2)``.
    ``Lasso`` needs ``max_iter=50000`` on this dataset -- the default 1000
    does not converge and emits a ``ConvergenceWarning``.
    """
    X_train, X_test, y_train, y_test = load_train_test()
    rows = []
    for alpha in alphas:
        lasso = Lasso(alpha=alpha, max_iter=50000).fit(X_train, y_train)
        ridge = Ridge(alpha=alpha).fit(X_train, y_train)
        lasso_zeros = int(np.sum(lasso.coef_ == 0))
        ridge_zeros = int(np.sum(ridge.coef_ == 0))
        lasso_r2 = round(r2_score(y_test, lasso.predict(X_test)), 4)
        ridge_r2 = round(r2_score(y_test, ridge.predict(X_test)), 4)
        rows.append((alpha, lasso_zeros, lasso_r2, ridge_zeros, ridge_r2))
    return rows


def lasso_cv_selection():
    """The alpha LassoCV picks by 5-fold cross-validation on the training
    split, how many coefficients it zeros, and which features it keeps.
    """
    X_train, X_test, y_train, y_test = load_train_test()
    model = LassoCV(cv=5, random_state=0, max_iter=50000).fit(X_train, y_train)
    zeros = int(np.sum(model.coef_ == 0))
    kept = [name for name, coef in zip(FEATURE_NAMES, model.coef_) if coef != 0]
    return {
        "alpha": float(model.alpha_),
        "zeros": zeros,
        "kept": kept,
        "r2": round(r2_score(y_test, model.predict(X_test)), 4),
    }


# --------------------------------------------------------------------------
# 2. The coefficient path: where each lasso coefficient hits exactly zero
# --------------------------------------------------------------------------


def coefficient_path(alphas):
    """Every coefficient, for both models, at every alpha in the sweep.

    Returns ``(lasso_path, ridge_path)``, each an array of shape
    ``(len(alphas), 10)`` in ``FEATURE_NAMES`` order, fitted on the full
    diabetes dataset (scaled) so the path is not split-dependent.
    """
    X, y = load_diabetes(return_X_y=True)
    lasso_path = np.zeros((len(alphas), X.shape[1]))
    ridge_path = np.zeros((len(alphas), X.shape[1]))
    for i, alpha in enumerate(alphas):
        lasso_path[i] = Lasso(alpha=alpha, max_iter=50000).fit(X, y).coef_
        ridge_path[i] = Ridge(alpha=alpha).fit(X, y).coef_
    return lasso_path, ridge_path


def alpha_where_each_lasso_coefficient_first_hits_zero(alphas):
    """The first alpha in the (ascending) sweep at which each lasso
    coefficient becomes exactly zero, and whether ridge ever hits zero at
    any alpha in the same sweep.

    Returns ``(zero_at, ridge_ever_zero)`` where ``zero_at`` maps each name
    in ``FEATURE_NAMES`` to the alpha value, or ``None`` if it never zeroed
    inside this sweep.
    """
    lasso_path, ridge_path = coefficient_path(alphas)
    zero_at = {name: None for name in FEATURE_NAMES}
    for i, alpha in enumerate(alphas):
        for j, name in enumerate(FEATURE_NAMES):
            if lasso_path[i, j] == 0.0 and zero_at[name] is None:
                zero_at[name] = float(alpha)
    ridge_ever_zero = bool(np.any(ridge_path == 0.0))
    return zero_at, ridge_ever_zero


# --------------------------------------------------------------------------
# 3. Does lasso pick the RIGHT features? A known sparse ground truth
# --------------------------------------------------------------------------


def sparse_recovery(alpha: float, noise: float, seed: int = 0):
    """Precision and recall of lasso's selected set against a KNOWN
    informative set, on synthetic data built with ``make_regression``.

    20 features, 5 informative, 200 rows. Returns
    ``(precision, recall, n_selected)``.
    """
    X, y, coef = make_regression(
        n_samples=200, n_features=20, n_informative=5, noise=noise, coef=True, random_state=seed
    )
    true_set = set(np.nonzero(coef)[0].tolist())
    model = Lasso(alpha=alpha, max_iter=50000).fit(X, y)
    picked = set(np.nonzero(model.coef_)[0].tolist())
    if not picked:
        return 0.0, 0.0, 0
    true_positives = len(true_set & picked)
    precision = round(true_positives / len(picked), 4)
    recall = round(true_positives / len(true_set), 4)
    return precision, recall, len(picked)


def sparse_recovery_grid(alphas, noises, seed: int = 0):
    """``sparse_recovery`` over every combination of alpha and noise.

    Returns rows of ``(alpha, noise, precision, recall, n_selected)``.
    """
    rows = []
    for noise in noises:
        for alpha in alphas:
            precision, recall, n_selected = sparse_recovery(alpha, noise, seed=seed)
            rows.append((alpha, noise, precision, recall, n_selected))
    return rows


def sparse_recovery_across_seeds(alpha: float, noise: float, seeds=range(10)):
    """``sparse_recovery`` averaged over several dataset seeds, so the
    finding does not rest on one lucky draw. Returns
    ``(mean_precision, mean_recall)``, each rounded to 4 places.
    """
    precisions = []
    recalls = []
    for seed in seeds:
        precision, recall, _n_selected = sparse_recovery(alpha, noise, seed=seed)
        precisions.append(precision)
        recalls.append(recall)
    return round(float(np.mean(precisions)), 4), round(float(np.mean(recalls)), 4)


# --------------------------------------------------------------------------
# 4. Scale-dependence: the penalty lives in whatever units the coefficients
#    happen to be in
# --------------------------------------------------------------------------


def scale_dependence(alpha: float = 1.0):
    """Lasso's selected feature set on the SAME data in three different
    units, at the SAME alpha: raw measurement units, standardised to unit
    variance, and scikit-learn's own bundled convention (unit L2 norm).

    Returns a dict with one entry per convention: ``kept`` (the feature
    names with a nonzero coefficient) and ``n_kept``.
    """
    X_raw, y = load_diabetes(return_X_y=True, scaled=False)
    X_unit_norm, _ = load_diabetes(return_X_y=True)  # sklearn's own scaled=True
    X_unit_variance = StandardScaler().fit_transform(X_raw)

    results = {}
    for label, X in (
        ("raw", X_raw),
        ("standardized", X_unit_variance),
        ("sklearn_unit_norm", X_unit_norm),
    ):
        model = Lasso(alpha=alpha, max_iter=50000).fit(X, y)
        kept = [name for name, coef in zip(FEATURE_NAMES, model.coef_) if coef != 0]
        results[label] = {"kept": kept, "n_kept": len(kept)}
    return results


# --------------------------------------------------------------------------
# 5. ElasticNet: the combination, and the alpha convention it does NOT share
#    with Ridge
# --------------------------------------------------------------------------


def elasticnet_sweep(alpha: float, l1_ratios):
    """ElasticNet's zero count and test R2 as l1_ratio moves from 0 (pure
    L2) to 1 (pure L1), at a fixed alpha, on the diabetes split.

    Returns rows of ``(l1_ratio, zeros, r2)``.
    """
    X_train, X_test, y_train, y_test = load_train_test()
    rows = []
    for l1_ratio in l1_ratios:
        model = ElasticNet(alpha=alpha, l1_ratio=l1_ratio, max_iter=50000).fit(X_train, y_train)
        zeros = int(np.sum(model.coef_ == 0))
        r2 = round(r2_score(y_test, model.predict(X_test)), 4)
        rows.append((l1_ratio, zeros, r2))
    return rows


def ridge_elasticnet_equivalence(alpha: float = 0.1):
    """Ridge and ElasticNet(l1_ratio=1) both define an L2-only penalty, but
    at DIFFERENT alpha scales: Ridge's objective sums the squared error,
    ElasticNet's averages it over n_samples. This measures the correction
    factor and confirms the two models agree once it is applied.

    Returns ``(ridge_coef_head, elasticnet_coef_head, max_abs_difference)``
    for the first three coefficients, after refitting Ridge at
    ``alpha * n_train``.
    """
    X_train, X_test, y_train, y_test = load_train_test()
    n_train = X_train.shape[0]
    ridge = Ridge(alpha=alpha * n_train).fit(X_train, y_train)
    elastic = ElasticNet(alpha=alpha, l1_ratio=0.0, max_iter=50000).fit(X_train, y_train)
    max_diff = float(np.max(np.abs(ridge.coef_ - elastic.coef_)))
    return ridge.coef_[:3].round(4).tolist(), elastic.coef_[:3].round(4).tolist(), round(max_diff, 4)


# --------------------------------------------------------------------------
# 6. Correlated predictors: ridge splits the weight, lasso picks one
# --------------------------------------------------------------------------


def near_duplicate_dataset(n: int = 300, seed: int = 0):
    """Two columns built from the same underlying signal plus tiny
    independent noise, so they are correlated at essentially 1.0, and a
    third column that is genuinely independent.
    """
    rng = np.random.default_rng(seed)
    base = rng.normal(size=n)
    x1 = base + rng.normal(scale=0.01, size=n)
    x2 = base + rng.normal(scale=0.01, size=n)
    x3 = rng.normal(size=n)
    X = np.column_stack([x1, x2, x3])
    y = 3.0 * x1 + 3.0 * x2 + 1.0 * x3 + rng.normal(scale=0.5, size=n)
    correlation = float(np.corrcoef(x1, x2)[0, 1])
    return X, y, correlation


def ridge_vs_lasso_on_duplicates(alphas):
    """Ridge and lasso coefficients on the near-duplicate dataset, at each
    alpha. Returns rows of ``(alpha, ridge_coefs, lasso_coefs)``, each a
    3-element list rounded to 4 places, in ``(x1, x2, x3)`` order.
    """
    X, y, _correlation = near_duplicate_dataset()
    rows = []
    for alpha in alphas:
        ridge = Ridge(alpha=alpha).fit(X, y)
        lasso = Lasso(alpha=alpha, max_iter=50000).fit(X, y)
        rows.append((alpha, ridge.coef_.round(4).tolist(), lasso.coef_.round(4).tolist()))
    return rows


# --------------------------------------------------------------------------
# 7. Ridge has a closed form; lasso does not
# --------------------------------------------------------------------------


def lasso_iteration_counts(alphas):
    """How many coordinate-descent iterations ``Lasso`` needed to converge,
    at each alpha, on the diabetes training split. Ridge has no equivalent
    attribute: it is solved directly by a single linear-algebra call, never
    iterated.
    """
    X_train, _X_test, y_train, _y_test = load_train_test()
    counts = {}
    for alpha in alphas:
        model = Lasso(alpha=alpha, max_iter=50000).fit(X_train, y_train)
        counts[alpha] = int(model.n_iter_)
    return counts


def ridge_has_no_iteration_count():
    """Confirm a fitted Ridge model carries no ``n_iter_`` under its default
    (closed-form) solver, unlike a fitted Lasso model.
    """
    X_train, _X_test, y_train, _y_test = load_train_test()
    ridge = Ridge(alpha=1.0).fit(X_train, y_train)
    lasso = Lasso(alpha=1.0, max_iter=50000).fit(X_train, y_train)
    return {
        "ridge_solver": ridge.solver,
        "ridge_has_n_iter": hasattr(ridge, "n_iter_") and ridge.n_iter_ is not None,
        "lasso_has_n_iter": hasattr(lasso, "n_iter_") and lasso.n_iter_ is not None,
    }


# --------------------------------------------------------------------------
# 8. The corner, in the smallest case that shows it: two correlated features
# --------------------------------------------------------------------------


def two_feature_dataset(n: int = 200, seed: int = 3):
    """Two strongly correlated features and a real, equal-weighted signal,
    small enough to draw. Returns ``(X, y, correlation)``.
    """
    rng = np.random.default_rng(seed)
    x1 = rng.normal(size=n)
    x2 = 0.9 * x1 + rng.normal(scale=0.3, size=n)
    X = np.column_stack([x1, x2])
    y = 2.0 * x1 + 2.0 * x2 + rng.normal(scale=1.0, size=n)
    correlation = float(np.corrcoef(x1, x2)[0, 1])
    return X, y, correlation


def two_feature_corner_demo(alphas):
    """Ridge and lasso coefficients on the two-feature dataset, at each
    alpha, alongside the unregularised (OLS) coefficients for reference.

    Returns ``(ols_coef, rows)`` where each row is
    ``(alpha, ridge_coefs, lasso_coefs)``, coefficients rounded to 4 places.
    """
    from sklearn.linear_model import LinearRegression

    X, y, _correlation = two_feature_dataset()
    ols = LinearRegression().fit(X, y).coef_.round(4).tolist()
    rows = []
    for alpha in alphas:
        ridge = Ridge(alpha=alpha).fit(X, y)
        lasso = Lasso(alpha=alpha, max_iter=50000).fit(X, y)
        rows.append((alpha, ridge.coef_.round(4).tolist(), lasso.coef_.round(4).tolist()))
    return ols, rows

#!/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 regularization_lib as r  # noqa: E402

ALPHA_GRID = [0.001, 0.01, 0.1, 1.0]
PATH_ALPHAS = np.logspace(-3, 2, 60)


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


def main() -> None:
    print("Day 151 -- regularization: ridge and lasso, measured")
    print("=" * 54)

    rule("1. Ridge never zeros; lasso zeros progressively more")
    print("   alpha   lasso-zeros   lasso-R2   ridge-zeros   ridge-R2")
    for alpha, lz, lr2, rz, rr2 in r.zero_counts_and_r2(ALPHA_GRID):
        print(f"  {alpha:6.3f}       {lz}/10      {lr2:.4f}        {rz}/10      {rr2:.4f}")

    rule("1b. LassoCV picks its own alpha")
    cv = r.lasso_cv_selection()
    print(f"  alpha={cv['alpha']:.5f}  zeros={cv['zeros']}/10  test R2={cv['r2']:.4f}")
    print(f"  kept: {cv['kept']}")

    rule("2. The coefficient path: where each lasso coefficient hits zero")
    zero_at, ridge_ever_zero = r.alpha_where_each_lasso_coefficient_first_hits_zero(PATH_ALPHAS)
    for name in r.FEATURE_NAMES:
        print(f"  {name:4s} zeros at alpha = {zero_at[name]:.4f}")
    print(f"  ridge ever zero across the same 60-point sweep: {ridge_ever_zero}")

    rule("3. Does lasso recover the RIGHT features? A known sparse truth")
    p1, r1, n1 = r.sparse_recovery(alpha=1.0, noise=1.0, seed=0)
    print(f"  alpha=1.0  noise=1.0   precision={p1:.4f}  recall={r1:.4f}  n_selected={n1}")
    p2, r2v, n2 = r.sparse_recovery(alpha=80.0, noise=30.0, seed=0)
    print(f"  alpha=80.0 noise=30.0  precision={p2:.4f}  recall={r2v:.4f}  n_selected={n2}")
    p3, r3, n3 = r.sparse_recovery(alpha=80.0, noise=10.0, seed=0)
    print(f"  alpha=80.0 noise=10.0  precision={p3:.4f}  recall={r3:.4f}  n_selected={n3}")
    mp_low, mr_low = r.sparse_recovery_across_seeds(alpha=1.0, noise=1.0)
    mp_high, mr_high = r.sparse_recovery_across_seeds(alpha=1.0, noise=10.0)
    print(f"  mean over 10 seeds, alpha=1.0 noise=1.0  : precision={mp_low:.4f}  recall={mr_low:.4f}")
    print(f"  mean over 10 seeds, alpha=1.0 noise=10.0 : precision={mp_high:.4f}  recall={mr_high:.4f}")

    rule("4. Regularization requires scaled features")
    scale = r.scale_dependence(alpha=1.0)
    for label in ("raw", "standardized", "sklearn_unit_norm"):
        print(f"  {label:18s} n_kept={scale[label]['n_kept']:2d}  kept={scale[label]['kept']}")

    rule("5. ElasticNet interpolates between ridge and lasso")
    print("  l1_ratio   zeros   R2")
    for l1_ratio, zeros, r2v in r.elasticnet_sweep(alpha=0.1, l1_ratios=[0.0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.0]):
        print(f"  {l1_ratio:6.1f}      {zeros}/10   {r2v:.4f}")

    rule("5b. Ridge and ElasticNet do not share an alpha scale")
    ridge_head, elastic_head, max_diff = r.ridge_elasticnet_equivalence(alpha=0.1)
    print(f"  Ridge(alpha=0.1 * n_train) coefs[:3]        : {ridge_head}")
    print(f"  ElasticNet(alpha=0.1, l1_ratio=0) coefs[:3] : {elastic_head}")
    print(f"  max abs difference after the n_train correction: {max_diff:.4f}")

    rule("6. Correlated predictors: ridge splits, lasso picks one")
    _X, _y, correlation = r.near_duplicate_dataset()
    print(f"  correlation between the two near-duplicate columns: {correlation:.6f}")
    for alpha, ridge_coefs, lasso_coefs in r.ridge_vs_lasso_on_duplicates([0.001, 0.1, 1.0, 10.0]):
        print(f"  alpha={alpha:6.3f}  ridge={ridge_coefs}  lasso={lasso_coefs}")

    rule("7. Ridge has a closed form; lasso needs iterations")
    info = r.ridge_has_no_iteration_count()
    print(f"  ridge solver: {info['ridge_solver']}  ridge has n_iter_: {info['ridge_has_n_iter']}")
    print(f"  lasso has n_iter_: {info['lasso_has_n_iter']}")
    counts = r.lasso_iteration_counts(ALPHA_GRID)
    for alpha, n_iter in counts.items():
        print(f"  lasso alpha={alpha:<6} n_iter_={n_iter}")

    rule("8. The corner: two correlated features, small enough to see it")
    ols, rows = r.two_feature_corner_demo([0.001, 0.5, 1.0, 3.0, 8.0])
    print(f"  OLS coefficients (no penalty): {ols}")
    for alpha, ridge_coefs, lasso_coefs in rows:
        print(f"  alpha={alpha:6.3f}  ridge={ridge_coefs}  lasso={lasso_coefs}")


if __name__ == "__main__":
    main()

"""Fourteen exercises in what a penalty actually does, measured.

Read `00_brief.md` first. Each function below is solved here and is a
`pytest.skip` in `starter/` naming exactly what to build and what to
assert. `regularization_lib.py` is complete -- it is the machinery, not the
exercise.

Run this suite on its own:

    .venv/bin/pytest examples -q

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

import numpy as np

import regularization_lib as r

ALPHA_GRID = [0.001, 0.01, 0.1, 1.0]
PATH_ALPHAS = np.logspace(-3, 2, 60)


def test_01_ridge_never_zeros_lasso_progressively_zeros():
    rows = r.zero_counts_and_r2(ALPHA_GRID)
    assert rows == [
        (0.001, 0, 0.3588, 0, 0.3586),
        (0.01, 1, 0.3541, 0, 0.3567),
        (0.1, 3, 0.355, 0, 0.369),
        (1.0, 8, 0.2782, 0, 0.357),
    ]
    ridge_zeros = [row[3] for row in rows]
    lasso_zeros = [row[1] for row in rows]
    assert ridge_zeros == [0, 0, 0, 0]
    # lasso's zero count never goes down as the penalty grows
    assert all(a <= b for a, b in zip(lasso_zeros, lasso_zeros[1:]))
    assert lasso_zeros[-1] > lasso_zeros[0]


def test_01b_lasso_cv_picks_alpha_and_six_features():
    result = r.lasso_cv_selection()
    assert round(result["alpha"], 5) == 0.07874
    assert result["zeros"] == 4
    assert result["kept"] == ["sex", "bmi", "bp", "s1", "s3", "s5"]
    assert result["r2"] == 0.3562


def test_02_the_coefficient_path_lasso_hits_exact_zeros_ridge_never_does():
    zero_at, ridge_ever_zero = r.alpha_where_each_lasso_coefficient_first_hits_zero(PATH_ALPHAS)
    # every lasso coefficient hits exactly zero somewhere in this sweep
    assert all(value is not None for value in zero_at.values())
    assert ridge_ever_zero is False


def test_02b_the_weakest_coefficients_zero_out_first():
    zero_at, _ridge_ever_zero = r.alpha_where_each_lasso_coefficient_first_hits_zero(PATH_ALPHAS)
    # s3 is the first coefficient lasso can afford to drop; bmi is the last
    weakest_first = min(zero_at, key=zero_at.get)
    strongest_last = max(zero_at, key=zero_at.get)
    assert weakest_first == "s3"
    assert strongest_last == "bmi"
    assert round(zero_at["s3"], 4) == 0.0032
    assert round(zero_at["bmi"], 4) == 2.4538


def test_03_lasso_recovers_the_right_features_at_low_noise():
    precision, recall, n_selected = r.sparse_recovery(alpha=1.0, noise=1.0, seed=0)
    assert precision == 1.0
    assert recall == 1.0
    assert n_selected == 5
    precision2, recall2, n_selected2 = r.sparse_recovery(alpha=0.1, noise=0.1, seed=0)
    assert precision2 == 1.0
    assert recall2 == 1.0
    assert n_selected2 == 5


def test_03b_recovery_degrades_with_noise_and_too_much_penalty():
    # honest failure: at high noise AND a heavy penalty, lasso can miss
    # most of the true informative set, or all of it
    precision, recall, n_selected = r.sparse_recovery(alpha=80.0, noise=30.0, seed=0)
    assert recall == 0.2
    assert n_selected == 1
    precision0, recall0, n_selected0 = r.sparse_recovery(alpha=80.0, noise=10.0, seed=0)
    assert recall0 == 0.0
    assert n_selected0 == 0
    # and this is not one unlucky seed: averaged over ten dataset seeds,
    # recovery is excellent at moderate noise and degrades at high noise
    mean_precision_low, mean_recall_low = r.sparse_recovery_across_seeds(alpha=1.0, noise=1.0)
    assert mean_precision_low == 1.0
    assert mean_recall_low == 0.98
    mean_precision_high, mean_recall_high = r.sparse_recovery_across_seeds(alpha=1.0, noise=10.0)
    assert mean_precision_high < mean_precision_low
    assert mean_precision_high == 0.6792


def test_04_regularization_requires_scaled_features():
    result = r.scale_dependence(alpha=1.0)
    assert result["raw"]["kept"] == [
        "age", "sex", "bmi", "bp", "s1", "s2", "s3", "s4", "s5", "s6",
    ]
    assert result["raw"]["n_kept"] == 10
    assert result["standardized"]["kept"] == ["sex", "bmi", "bp", "s1", "s3", "s5", "s6"]
    assert result["standardized"]["n_kept"] == 7
    assert result["sklearn_unit_norm"]["kept"] == ["bmi", "bp", "s5"]
    assert result["sklearn_unit_norm"]["n_kept"] == 3
    # three different answers, same data, same alpha -- only the units differ
    counts = {result[key]["n_kept"] for key in result}
    assert len(counts) == 3


def test_05_elasticnet_interpolates_between_ridge_and_lasso():
    rows = r.elasticnet_sweep(alpha=0.1, l1_ratios=[0.0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.0])
    assert rows == [
        (0.0, 0, 0.0555),
        (0.1, 0, 0.0605),
        (0.3, 0, 0.0741),
        (0.5, 0, 0.0963),
        (0.7, 1, 0.1389),
        (0.9, 0, 0.2511),
        (1.0, 3, 0.355),
    ]
    # pure L1 (l1_ratio=1.0) is the only setting that matches plain Lasso's
    # own zero count at the same alpha
    lasso_only = r.zero_counts_and_r2([0.1])[0]
    assert rows[-1][1] == lasso_only[1]


def test_05b_ridge_and_elasticnet_do_not_share_an_alpha_scale():
    ridge_head, elastic_head, max_diff = r.ridge_elasticnet_equivalence(alpha=0.1)
    # naively comparing Ridge(alpha=0.1) to ElasticNet(alpha=0.1, l1_ratio=0)
    # would NOT agree -- Ridge sums the squared error, ElasticNet averages
    # it, so the alphas differ by a factor of n_train. Once that correction
    # is applied (Ridge fitted at alpha * n_train), the two match closely.
    assert max_diff < 0.001
    assert ridge_head == [6.3098, 1.229, 22.6076]
    assert elastic_head == [6.3097, 1.229, 22.6076]


def test_06_ridge_splits_the_weight_between_near_duplicates():
    _X, _y, correlation = r.near_duplicate_dataset()
    assert correlation > 0.999
    rows = r.ridge_vs_lasso_on_duplicates([1.0])
    _alpha, ridge_coefs, lasso_coefs = rows[0]
    ridge_x1, ridge_x2, _ridge_x3 = ridge_coefs
    lasso_x1, lasso_x2, _lasso_x3 = lasso_coefs
    # ridge splits the true combined weight of 6.0 roughly evenly
    assert abs(ridge_x1 - ridge_x2) < 0.15
    assert abs((ridge_x1 + ridge_x2) - 6.0) < 0.2
    # lasso picks one and drives the other to exactly zero
    assert lasso_x2 == 0.0
    assert lasso_x1 > 5.0


def test_06b_enough_penalty_zeros_both_duplicates_in_lasso_but_ridge_still_splits():
    rows = r.ridge_vs_lasso_on_duplicates([10.0])
    _alpha, ridge_coefs, lasso_coefs = rows[0]
    ridge_x1, ridge_x2, _ridge_x3 = ridge_coefs
    assert lasso_coefs[0] == 0.0
    assert lasso_coefs[1] == 0.0
    # ridge, at the same alpha, still has both near-duplicate coefficients
    # alive and still close to each other
    assert ridge_x1 > 2.5
    assert ridge_x2 > 2.5
    assert abs(ridge_x1 - ridge_x2) < 0.15


def test_07_ridge_has_a_closed_form_lasso_needs_iterations():
    info = r.ridge_has_no_iteration_count()
    assert info["ridge_has_n_iter"] is False
    assert info["lasso_has_n_iter"] is True
    counts = r.lasso_iteration_counts(ALPHA_GRID)
    assert counts == {0.001: 368, 0.01: 62, 0.1: 135, 1.0: 6}
    # every count is real work, and none of them hit the ceiling
    assert all(0 < n < 50000 for n in counts.values())


def test_08_the_corner_two_correlated_features_and_a_lasso_zero():
    ols, rows = r.two_feature_corner_demo([0.001, 0.5, 1.0, 3.0, 8.0])
    assert ols == [1.9564, 1.9381]
    by_alpha = {alpha: (ridge, lasso) for alpha, ridge, lasso in rows}
    # at alpha=3.0, lasso's solution has landed exactly on the axis --
    # the second coefficient is exactly zero -- while ridge's has not
    ridge_at_3, lasso_at_3 = by_alpha[3.0]
    assert lasso_at_3 == [0.8919, 0.0]
    assert ridge_at_3[0] != 0.0 and ridge_at_3[1] != 0.0
    assert abs(ridge_at_3[1]) > 1.5


def test_08b_at_a_tiny_alpha_both_models_agree_with_ols():
    ols, rows = r.two_feature_corner_demo([0.001, 8.0])
    by_alpha = {alpha: (ridge, lasso) for alpha, ridge, lasso in rows}
    ridge_tiny, lasso_tiny = by_alpha[0.001]
    for ols_coef, ridge_coef in zip(ols, ridge_tiny):
        assert abs(ols_coef - ridge_coef) < 0.01
    for ols_coef, lasso_coef in zip(ols, lasso_tiny):
        assert abs(ols_coef - lasso_coef) < 0.01
    # and at a large enough alpha, lasso has zeroed everything while
    # ridge, which never zeros, has merely shrunk
    ridge_big, lasso_big = by_alpha[8.0]
    assert lasso_big == [0.0, 0.0]
    assert ridge_big[0] != 0.0 and ridge_big[1] != 0.0

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

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

import numpy as np

import regularization_lib as r


def test_the_diabetes_split_is_reproducible_and_shaped_right():
    a = r.load_train_test()
    b = r.load_train_test()
    X_train_a, X_test_a, y_train_a, y_test_a = a
    X_train_b, _X_test_b, _y_train_b, _y_test_b = b
    assert X_train_a.shape == (331, 10)
    assert X_test_a.shape == (111, 10)
    assert y_train_a.shape == (331,)
    assert y_test_a.shape == (111,)
    # same seed, same split -- every downstream measurement depends on this
    assert np.array_equal(X_train_a, X_train_b)


def test_the_near_duplicate_columns_really_are_near_duplicate():
    X, y, correlation = r.near_duplicate_dataset()
    assert X.shape == (300, 3)
    assert y.shape == (300,)
    assert correlation > 0.999
    # the third column is not part of the near-duplicate pair
    assert abs(float(np.corrcoef(X[:, 0], X[:, 2])[0, 1])) < 0.3


def test_the_synthetic_sparse_dataset_really_has_five_informative_features():
    from sklearn.datasets import make_regression

    X, y, coef = make_regression(
        n_samples=200, n_features=20, n_informative=5, noise=1.0, coef=True, random_state=0
    )
    assert X.shape == (200, 20)
    assert int(np.sum(coef != 0)) == 5


def test_ridge_solves_directly_and_lasso_iterates():
    info = r.ridge_has_no_iteration_count()
    assert info["ridge_has_n_iter"] is False
    assert info["lasso_has_n_iter"] is True
    counts = r.lasso_iteration_counts([0.001, 1.0])
    # every count is a positive number of iterations, and none of them
    # hit the ceiling -- so no ConvergenceWarning was silently swallowed
    for alpha, n_iter in counts.items():
        assert 0 < n_iter < 50000, f"alpha={alpha} iterated {n_iter} times"

lesson_id: D151
day: 151
kind: guided-build
languages:
  - python
  - bash
setup_commands:
  - cd labs/sections/machine-learning/day-151-regularization-ridge-and-lasso
  - python3 -m venv .venv
  - .venv/bin/pip install -r requirements/requirements.txt
  - >-
    .venv/bin/python3 -c "import numpy, sklearn; print(numpy.__version__,
    sklearn.__version__)"
run_commands:
  - .venv/bin/pytest examples -q
  - .venv/bin/pytest starter -q
  - .venv/bin/python3 examples/report_measurements.py
test_commands:
  - bash tests/run_tests.sh
cleanup_commands:
  - >-
    find . -path ./.venv -prune -o -type d -name '__pycache__' -print -exec rm -rf -- {}
    +
  - rm -rf .pytest_cache
  - 'rm -rf .venv  # optional: removes the lab virtual environment'
  - 'git checkout -- starter/  # optional: reset your work'
requires_network: true
requires_api_key: false
estimated_minutes: 60
last_executed: '2026-08-27'
executed_on: >-
  macOS 26.5.2 (Apple Silicon, arm64, CPU only -- no GPU is needed or used), Python
  3.14.0, numpy 2.5.2, scikit-learn 1.9.0, pytest 9.1.1, bash 3.2.57 -- bash
  tests/run_tests.sh -> 14 checks, 0 failure(s), exit 0. pytest examples -q -> 18
  passed. pytest starter -q -> 4 passed, 14 skipped (the four machinery checks in
  test_regularization_lib.py are solved in both directories; the fourteen exercise
  stubs in starter/test_regularization_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 -- load_diabetes ships bundled inside the scikit-learn package
  itself and is read from local disk, never fetched at runtime; every synthetic
  dataset is generated on the spot from a seeded numpy.random.default_rng or
  sklearn.datasets.make_regression; harness check 9 confirms no URL appears anywhere
  in starter/ or examples/ source. Section 7 of the harness copies examples/ into a
  mktemp-d scratch directory, confirms 18 passed, rewrites `assert lasso_x2 == 0.0` in
  test_06_ridge_splits_the_weight_between_near_duplicates to `assert lasso_x2 ==
  999.0`, confirms a non-zero exit naming the failing test, and removes the scratch
  directory. Separately, by hand, `assert result["raw"]["n_kept"] == 10` was changed
  to 3 in examples/test_regularization_claims.py and the whole harness re-run: it
  reported 14 checks, 2 failure(s) and exited 1 (both the pytest run and the
  pytest-free direct reproduction in section 2 caught it); the file was restored and
  the harness returned to 14 checks, 0 failure(s), exit 0. MEASURED PAIRS, all
  captured verbatim in expected-output/measured-values.txt. (1) HEADLINE, on
  load_diabetes(return_X_y=True) with train_test_split(test_size=0.25,
  random_state=0): at alpha 0.001, 0.01, 0.1 and 1.0, lasso zeros 0, 1, 3 and 8 of 10
  coefficients (test R2 0.3588, 0.3541, 0.3550, 0.2782) while ridge zeros 0 of 10 at
  every one of those alphas (test R2 0.3586, 0.3567, 0.3690, 0.3570). LassoCV(cv=5)
  picks alpha=0.07874, zeros 4/10, and keeps ['sex','bmi','bp','s1','s3','s5'] with
  test R2 0.3562. (2) THE COEFFICIENT PATH, swept over 60 log-spaced alphas from 0.001
  to 100 on the full dataset: every one of the ten lasso coefficients hits EXACTLY
  zero somewhere in the sweep -- s3 first, at alpha=0.0032, bmi last, at alpha=2.4538
  -- while not one ridge coefficient hits zero at any alpha in the same sweep. (3)
  SPARSE RECOVERY against a KNOWN ground truth (make_regression, 5 informative of 20
  features): at alpha=1.0, noise=1.0, lasso recovers precision=1.0, recall=1.0,
  selecting exactly the 5 true features. At alpha=80.0 (too much penalty) with
  noise=30.0, recall falls to 0.2 (1 of 5 recovered); at noise=10.0, recall falls to
  0.0 (the penalty zeroed everything). Averaged over 10 dataset seeds at alpha=1.0:
  mean precision 1.0 at noise=1.0, falling to mean precision 0.6792 at noise=10.0,
  both with mean recall 0.98. (4) SCALE-DEPENDENCE, THE MOST PRACTICALLY IMPORTANT
  RESULT: at the identical alpha=1.0, on the identical diabetes data, lasso in raw
  measurement units keeps all 10 of 10 features; standardised to unit variance
  (StandardScaler) it keeps 7 (['sex','bmi','bp','s1','s3','s5','s6']); in
  scikit-learn's own bundled "scaled=True" convention (unit L2 norm, not unit
  variance) it keeps 3 (['bmi','bp','s5']). Three different answers, same data, same
  alpha, same random_state -- because a lasso penalty is applied in whatever units the
  coefficients happen to be in, and even "scaled" is not one convention. (5)
  ELASTICNET at alpha=0.1: zero count and test R2 move from (0 zeros, R2=0.0555) at
  l1_ratio=0.0 up to (3 zeros, R2=0.3550) at l1_ratio=1.0, and the l1_ratio=1.0 row's
  zero count and R2 match plain Lasso(alpha=0.1) exactly. (6) RIDGE AND ELASTICNET DO
  NOT SHARE AN ALPHA SCALE: Ridge's objective sums the squared error over all rows;
  ElasticNet's averages it over n_samples. Naively comparing Ridge(alpha=0.1) to
  ElasticNet(alpha=0.1, l1_ratio=0.0) would NOT agree; fitting Ridge at
  alpha*n_train=33.1 instead matches ElasticNet(alpha=0.1, l1_ratio=0.0) to within
  0.0001 on the first three coefficients (measured max abs difference 3.8e-6). (7)
  CORRELATED PREDICTORS (two columns built from one shared signal, correlation
  0.999918): at alpha=1.0, ridge splits the true combined weight of 6.0 almost evenly
  (3.048 and 2.9742), while lasso keeps one at 5.0848 and sets the other to EXACTLY
  0.0. At alpha=10.0, lasso has zeroed BOTH
  duplicate coefficients while ridge STILL splits its two nonzero coefficients almost
  evenly (2.9724 and 2.9646) -- ridge never produces the corner, at any alpha tried.
  (8) RIDGE HAS A CLOSED FORM; LASSO DOES NOT: a fitted Ridge model carries no
  n_iter_ attribute at all (solved by one linear-algebra call, solver="auto");
  a fitted Lasso model does, and the iteration count varies with alpha: 368 at
  alpha=0.001, 62 at alpha=0.01, 135 at alpha=0.1, 6 at alpha=1.0 -- none of them
  hitting the max_iter=50000 ceiling, so no ConvergenceWarning was silently
  swallowed. (9) THE CORNER, in the smallest case that shows it (two features
  correlated at 0.9999, OLS solution [1.9564, 1.9381]): as alpha rises through 0.001,
  0.5, 1.0, 3.0 and 8.0, ridge's two coefficients shrink together and stay both
  nonzero the whole way (1.9345 and 1.8820 even at alpha=8.0); lasso's second
  coefficient lands on EXACTLY 0.0 at alpha=3.0, while its first coefficient is still
  0.8919 -- the corner of the L1 diamond sitting on an axis. FOUR HONESTY CALLS.
  FIRST: exercise 3b shows lasso can fail outright -- at high noise and a heavy
  penalty, recall drops to 0.0, meaning the penalty zeroed every informative feature.
  Regularization is not free precision; too much of it destroys the signal it was
  supposed to isolate. SECOND: the near-duplicate columns in exercise 6 are
  correlated at 0.9999, not identically collinear, so lasso's split is not guaranteed
  to be an EXACT hard zero at every dataset seed -- harness check 8 re-runs the
  construction at three more seeds and found seed 3 gives lasso coefficients
  [5.0215, 0.0007] rather than an exact [x, 0.0]; the assertion there was loosened to
  "a 20x-or-greater asymmetric split" rather than an exact zero, because that is what
  the measurement actually supports beyond the one quoted seed. THIRD: the
  scale-dependence result (exercise 4) does not merely differ in HOW MANY features
  survive -- StandardScaler's unit-variance convention and scikit-learn's own
  unit-L2-norm convention are BOTH legitimately called "scaled" and give different
  answers (7 vs 3 features kept) at the identical nominal alpha, which means alpha
  values are not portable even between two reasonable definitions of "scaled." FOURTH:
  the Ridge/ElasticNet alpha-scale mismatch (exercise 5b) is a real API gotcha, not a
  contrived one -- an engineer who naively swaps ElasticNet(l1_ratio=0) in for Ridge
  without correcting alpha by n_train will silently regularize far more aggressively
  than intended, which is exactly what the uncorrected comparison in this lab's own
  exploratory work initially produced (R2 0.0555 against Ridge's 0.3690 at a
  nominally-matched alpha=0.1).

# 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 — scipy 1.18.1, joblib
1.5.3, threadpoolctl 3.6.0 — are recorded in `../expected-output/FIELDS.md`.

## Why the versions are pinned exactly

Every number in this lab comes from `sklearn.linear_model.Ridge`,
`Lasso`, `LassoCV` and `ElasticNet` fitted on `sklearn.datasets.load_diabetes`
and `sklearn.datasets.make_regression`. Coordinate-descent convergence
paths, default solver choices, and exact floating-point results can shift
between scikit-learn releases even when the documented behaviour does not.

What does not depend on the pins: the directions of every result — ridge
never zeroing a coefficient at any alpha tried, lasso zeroing progressively
more, the constraint-region geometry that makes a lasso coefficient land
exactly on zero, the alpha-scale mismatch between `Ridge` and `ElasticNet`,
and the fact that a penalty applied to unscaled features selects a
different set than the same penalty applied to scaled ones.
`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
only dataset downloaded from anywhere is `load_diabetes`, which ships
bundled inside the scikit-learn package itself and is read from local
disk, never fetched over the network at runtime. The synthetic datasets
are generated on the spot with a seeded generator.

## 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 151 lab brief — What the Penalty Does

Day 145 already measured that a ridge penalty rescued an overfit
degree-24 polynomial by a factor of 39,588, and that a penalty raises
training error while lowering the variance of what the model learns.

That day treated "regularization" as one thing. It is not. This lab
measures the contrast between the two most common penalties — L2 (ridge)
and L1 (lasso) — on the same dataset, with the same alphas, so the
difference is not a claim you take on faith.

## The claim you are here to measure

> Ridge shrinks. Lasso shrinks and selects. The difference is not a
> matter of degree — it is a difference in shape, and the shape has a
> geometric reason.

Exercise 1 puts the two side by side on `sklearn.datasets.load_diabetes`,
sweeping alpha from 0.001 to 1.0:

| alpha | lasso zeros | lasso R2 | ridge zeros | ridge R2 |
| --- | --- | --- | --- | --- |
| 0.001 | 0/10 | 0.3588 | 0/10 | 0.3586 |
| 0.01 | 1/10 | 0.3541 | 0/10 | 0.3567 |
| 0.1 | 3/10 | 0.3550 | 0/10 | 0.3690 |
| 1.0 | 8/10 | 0.2782 | 0/10 | 0.3570 |

Read the ridge column first. **It never moves.** Zero coefficients at
every alpha tried, from barely-there to aggressive. Ridge shrinks
everything toward zero and never quite arrives.

Now read the lasso column. It climbs steadily to 8 of 10 coefficients
zeroed. Lasso is doing feature selection, not just shrinkage — and it is
doing it as a side effect of the penalty shape, not because anyone told
it which features matter.

## Why: the geometry, made numeric

Exercise 8 builds the smallest case that shows the mechanism: two
strongly correlated features (correlation 0.9999) and an ordinary
least-squares solution of `[1.9564, 1.9381]`. As alpha grows:

| alpha | ridge | lasso |
| --- | --- | --- |
| 0.001 | `[1.9564, 1.9381]` | `[1.9583, 1.9352]` |
| 3.0 | `[1.9507, 1.9141]` | `[0.8919, 0.0]` |
| 8.0 | `[1.9345, 1.882]` | `[0.0, 0.0]` |

At alpha=3.0, lasso's second coefficient is **exactly** 0.0 — not small,
not rounded, exactly zero — while ridge's is still 1.91. That is the
corner of the L1 constraint region landing on an axis. A diamond has
corners on its axes; a circle (ridge's constraint region) does not
touch an axis except at a single tangent point that requires infinite
penalty to reach. Sections "How it works" and the architecture diagram
in the lesson walk through why the shape of the constraint region
forces this.

## What that buys you, and what it costs you

**Exercise 3** confirms lasso does not merely shrink toward *some*
sparse answer — on `make_regression` with a known set of 5 informative
features out of 20, it recovers precisely that set (precision 1.0,
recall 1.0) at a sensible alpha. **Exercise 3b** is the honest half:
push the alpha too high on noisy data and lasso can zero out the truth
entirely (recall drops to 0.0). More penalty is not free precision.

**Exercise 4** is the one with the most practical consequence. The same
alpha, on the same data, in three different units — raw measurement
units, standardized to unit variance, and scikit-learn's own bundled
convention (unit L2 norm) — selects **10, 7, and 3 features
respectively**. The penalty lives in whatever units the coefficients
happen to be in. Skip the scaling step and you have not actually
regularized anything meaningful; you have penalized whichever feature
happened to have small natural units.

**Exercise 6** connects straight back to Day 150's multicollinearity: on
two near-duplicate columns, ridge splits the combined weight roughly
evenly between them (each gets about half); lasso picks one and zeros
the other, arbitrarily as far as the data is concerned.

**Exercise 7** is the practical cost of the shape: ridge is one
linear-algebra call (no `n_iter_` attribute at all); lasso is solved
iteratively — coordinate descent — because its penalty is not
differentiable at zero and no closed form exists.

## How to work

1. Build the environment (see the lab `README.md`).
2. Run `.venv/bin/pytest starter -q`. You will see four passes (the
   machinery checks in `test_regularization_lib.py`) and fourteen skips.
3. Replace one `pytest.skip(...)` at a time with real code. The skip
   text names the exact helper and the exact value to assert.
4. Print every measured pair. 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 `regularization_lib.py`, `test_regularization_lib.py` and
`test_regularization_claims.py`; pytest aborts on the module-name
collision. Run them separately, always.

## The scope of this lab

This lab does not re-teach that regularization trades variance for
bias — Day 145 measured that. It does not build linear regression from
scratch — Day 148 and Day 149 own that. It does not cover which metric
to report — Day 152 owns that. It measures exactly one thing in depth:
what the L1 and L2 penalty shapes do differently, and why.

"""Ridge and lasso, measured: what the penalty does, and why the two differ
in kind rather than degree.

Day 145 already measured that a ridge penalty rescued an overfit polynomial
by a factor of 39,588, and that a penalty trades variance for bias. This
module does not re-measure that trade; it measures the CONTRAST between the
two penalty shapes: an L2 penalty that shrinks every coefficient toward zero
but never quite reaches it, and an L1 penalty that drives some coefficients
to exactly zero and leaves the rest alone.

Everything here is deterministic given a seed, and every dataset is either
scikit-learn's bundled ``load_diabetes`` or generated on the spot with
``numpy.random.default_rng`` or ``sklearn.datasets.make_regression``.
"""

from __future__ import annotations

import numpy as np

from sklearn.datasets import load_diabetes, make_regression
from sklearn.linear_model import ElasticNet, Lasso, LassoCV, Ridge
from sklearn.metrics import r2_score
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

FEATURE_NAMES = ["age", "sex", "bmi", "bp", "s1", "s2", "s3", "s4", "s5", "s6"]


# --------------------------------------------------------------------------
# 1. Ridge zeros nothing; lasso zeros progressively more
# --------------------------------------------------------------------------


def load_train_test(scaled: bool = True, seed: int = 0):
    """The diabetes dataset, split the way every exercise in this lab uses.

    ``scaled=True`` is scikit-learn's own bundled version: mean-centred and
    scaled so every column has unit L2 norm. ``scaled=False`` returns the raw
    measurement units (age in years, bmi as a ratio, bp in mm Hg, and so on).
    """
    X, y = load_diabetes(return_X_y=True, scaled=scaled)
    return train_test_split(X, y, test_size=0.25, random_state=seed)


def zero_counts_and_r2(alphas):
    """For each alpha, how many lasso coefficients are exactly zero, how
    many ridge coefficients are exactly zero, and each model's test R2.

    Returns rows of ``(alpha, lasso_zeros, lasso_r2, ridge_zeros, ridge_r2)``.
    ``Lasso`` needs ``max_iter=50000`` on this dataset -- the default 1000
    does not converge and emits a ``ConvergenceWarning``.
    """
    X_train, X_test, y_train, y_test = load_train_test()
    rows = []
    for alpha in alphas:
        lasso = Lasso(alpha=alpha, max_iter=50000).fit(X_train, y_train)
        ridge = Ridge(alpha=alpha).fit(X_train, y_train)
        lasso_zeros = int(np.sum(lasso.coef_ == 0))
        ridge_zeros = int(np.sum(ridge.coef_ == 0))
        lasso_r2 = round(r2_score(y_test, lasso.predict(X_test)), 4)
        ridge_r2 = round(r2_score(y_test, ridge.predict(X_test)), 4)
        rows.append((alpha, lasso_zeros, lasso_r2, ridge_zeros, ridge_r2))
    return rows


def lasso_cv_selection():
    """The alpha LassoCV picks by 5-fold cross-validation on the training
    split, how many coefficients it zeros, and which features it keeps.
    """
    X_train, X_test, y_train, y_test = load_train_test()
    model = LassoCV(cv=5, random_state=0, max_iter=50000).fit(X_train, y_train)
    zeros = int(np.sum(model.coef_ == 0))
    kept = [name for name, coef in zip(FEATURE_NAMES, model.coef_) if coef != 0]
    return {
        "alpha": float(model.alpha_),
        "zeros": zeros,
        "kept": kept,
        "r2": round(r2_score(y_test, model.predict(X_test)), 4),
    }


# --------------------------------------------------------------------------
# 2. The coefficient path: where each lasso coefficient hits exactly zero
# --------------------------------------------------------------------------


def coefficient_path(alphas):
    """Every coefficient, for both models, at every alpha in the sweep.

    Returns ``(lasso_path, ridge_path)``, each an array of shape
    ``(len(alphas), 10)`` in ``FEATURE_NAMES`` order, fitted on the full
    diabetes dataset (scaled) so the path is not split-dependent.
    """
    X, y = load_diabetes(return_X_y=True)
    lasso_path = np.zeros((len(alphas), X.shape[1]))
    ridge_path = np.zeros((len(alphas), X.shape[1]))
    for i, alpha in enumerate(alphas):
        lasso_path[i] = Lasso(alpha=alpha, max_iter=50000).fit(X, y).coef_
        ridge_path[i] = Ridge(alpha=alpha).fit(X, y).coef_
    return lasso_path, ridge_path


def alpha_where_each_lasso_coefficient_first_hits_zero(alphas):
    """The first alpha in the (ascending) sweep at which each lasso
    coefficient becomes exactly zero, and whether ridge ever hits zero at
    any alpha in the same sweep.

    Returns ``(zero_at, ridge_ever_zero)`` where ``zero_at`` maps each name
    in ``FEATURE_NAMES`` to the alpha value, or ``None`` if it never zeroed
    inside this sweep.
    """
    lasso_path, ridge_path = coefficient_path(alphas)
    zero_at = {name: None for name in FEATURE_NAMES}
    for i, alpha in enumerate(alphas):
        for j, name in enumerate(FEATURE_NAMES):
            if lasso_path[i, j] == 0.0 and zero_at[name] is None:
                zero_at[name] = float(alpha)
    ridge_ever_zero = bool(np.any(ridge_path == 0.0))
    return zero_at, ridge_ever_zero


# --------------------------------------------------------------------------
# 3. Does lasso pick the RIGHT features? A known sparse ground truth
# --------------------------------------------------------------------------


def sparse_recovery(alpha: float, noise: float, seed: int = 0):
    """Precision and recall of lasso's selected set against a KNOWN
    informative set, on synthetic data built with ``make_regression``.

    20 features, 5 informative, 200 rows. Returns
    ``(precision, recall, n_selected)``.
    """
    X, y, coef = make_regression(
        n_samples=200, n_features=20, n_informative=5, noise=noise, coef=True, random_state=seed
    )
    true_set = set(np.nonzero(coef)[0].tolist())
    model = Lasso(alpha=alpha, max_iter=50000).fit(X, y)
    picked = set(np.nonzero(model.coef_)[0].tolist())
    if not picked:
        return 0.0, 0.0, 0
    true_positives = len(true_set & picked)
    precision = round(true_positives / len(picked), 4)
    recall = round(true_positives / len(true_set), 4)
    return precision, recall, len(picked)


def sparse_recovery_grid(alphas, noises, seed: int = 0):
    """``sparse_recovery`` over every combination of alpha and noise.

    Returns rows of ``(alpha, noise, precision, recall, n_selected)``.
    """
    rows = []
    for noise in noises:
        for alpha in alphas:
            precision, recall, n_selected = sparse_recovery(alpha, noise, seed=seed)
            rows.append((alpha, noise, precision, recall, n_selected))
    return rows


def sparse_recovery_across_seeds(alpha: float, noise: float, seeds=range(10)):
    """``sparse_recovery`` averaged over several dataset seeds, so the
    finding does not rest on one lucky draw. Returns
    ``(mean_precision, mean_recall)``, each rounded to 4 places.
    """
    precisions = []
    recalls = []
    for seed in seeds:
        precision, recall, _n_selected = sparse_recovery(alpha, noise, seed=seed)
        precisions.append(precision)
        recalls.append(recall)
    return round(float(np.mean(precisions)), 4), round(float(np.mean(recalls)), 4)


# --------------------------------------------------------------------------
# 4. Scale-dependence: the penalty lives in whatever units the coefficients
#    happen to be in
# --------------------------------------------------------------------------


def scale_dependence(alpha: float = 1.0):
    """Lasso's selected feature set on the SAME data in three different
    units, at the SAME alpha: raw measurement units, standardised to unit
    variance, and scikit-learn's own bundled convention (unit L2 norm).

    Returns a dict with one entry per convention: ``kept`` (the feature
    names with a nonzero coefficient) and ``n_kept``.
    """
    X_raw, y = load_diabetes(return_X_y=True, scaled=False)
    X_unit_norm, _ = load_diabetes(return_X_y=True)  # sklearn's own scaled=True
    X_unit_variance = StandardScaler().fit_transform(X_raw)

    results = {}
    for label, X in (
        ("raw", X_raw),
        ("standardized", X_unit_variance),
        ("sklearn_unit_norm", X_unit_norm),
    ):
        model = Lasso(alpha=alpha, max_iter=50000).fit(X, y)
        kept = [name for name, coef in zip(FEATURE_NAMES, model.coef_) if coef != 0]
        results[label] = {"kept": kept, "n_kept": len(kept)}
    return results


# --------------------------------------------------------------------------
# 5. ElasticNet: the combination, and the alpha convention it does NOT share
#    with Ridge
# --------------------------------------------------------------------------


def elasticnet_sweep(alpha: float, l1_ratios):
    """ElasticNet's zero count and test R2 as l1_ratio moves from 0 (pure
    L2) to 1 (pure L1), at a fixed alpha, on the diabetes split.

    Returns rows of ``(l1_ratio, zeros, r2)``.
    """
    X_train, X_test, y_train, y_test = load_train_test()
    rows = []
    for l1_ratio in l1_ratios:
        model = ElasticNet(alpha=alpha, l1_ratio=l1_ratio, max_iter=50000).fit(X_train, y_train)
        zeros = int(np.sum(model.coef_ == 0))
        r2 = round(r2_score(y_test, model.predict(X_test)), 4)
        rows.append((l1_ratio, zeros, r2))
    return rows


def ridge_elasticnet_equivalence(alpha: float = 0.1):
    """Ridge and ElasticNet(l1_ratio=1) both define an L2-only penalty, but
    at DIFFERENT alpha scales: Ridge's objective sums the squared error,
    ElasticNet's averages it over n_samples. This measures the correction
    factor and confirms the two models agree once it is applied.

    Returns ``(ridge_coef_head, elasticnet_coef_head, max_abs_difference)``
    for the first three coefficients, after refitting Ridge at
    ``alpha * n_train``.
    """
    X_train, X_test, y_train, y_test = load_train_test()
    n_train = X_train.shape[0]
    ridge = Ridge(alpha=alpha * n_train).fit(X_train, y_train)
    elastic = ElasticNet(alpha=alpha, l1_ratio=0.0, max_iter=50000).fit(X_train, y_train)
    max_diff = float(np.max(np.abs(ridge.coef_ - elastic.coef_)))
    return ridge.coef_[:3].round(4).tolist(), elastic.coef_[:3].round(4).tolist(), round(max_diff, 4)


# --------------------------------------------------------------------------
# 6. Correlated predictors: ridge splits the weight, lasso picks one
# --------------------------------------------------------------------------


def near_duplicate_dataset(n: int = 300, seed: int = 0):
    """Two columns built from the same underlying signal plus tiny
    independent noise, so they are correlated at essentially 1.0, and a
    third column that is genuinely independent.
    """
    rng = np.random.default_rng(seed)
    base = rng.normal(size=n)
    x1 = base + rng.normal(scale=0.01, size=n)
    x2 = base + rng.normal(scale=0.01, size=n)
    x3 = rng.normal(size=n)
    X = np.column_stack([x1, x2, x3])
    y = 3.0 * x1 + 3.0 * x2 + 1.0 * x3 + rng.normal(scale=0.5, size=n)
    correlation = float(np.corrcoef(x1, x2)[0, 1])
    return X, y, correlation


def ridge_vs_lasso_on_duplicates(alphas):
    """Ridge and lasso coefficients on the near-duplicate dataset, at each
    alpha. Returns rows of ``(alpha, ridge_coefs, lasso_coefs)``, each a
    3-element list rounded to 4 places, in ``(x1, x2, x3)`` order.
    """
    X, y, _correlation = near_duplicate_dataset()
    rows = []
    for alpha in alphas:
        ridge = Ridge(alpha=alpha).fit(X, y)
        lasso = Lasso(alpha=alpha, max_iter=50000).fit(X, y)
        rows.append((alpha, ridge.coef_.round(4).tolist(), lasso.coef_.round(4).tolist()))
    return rows


# --------------------------------------------------------------------------
# 7. Ridge has a closed form; lasso does not
# --------------------------------------------------------------------------


def lasso_iteration_counts(alphas):
    """How many coordinate-descent iterations ``Lasso`` needed to converge,
    at each alpha, on the diabetes training split. Ridge has no equivalent
    attribute: it is solved directly by a single linear-algebra call, never
    iterated.
    """
    X_train, _X_test, y_train, _y_test = load_train_test()
    counts = {}
    for alpha in alphas:
        model = Lasso(alpha=alpha, max_iter=50000).fit(X_train, y_train)
        counts[alpha] = int(model.n_iter_)
    return counts


def ridge_has_no_iteration_count():
    """Confirm a fitted Ridge model carries no ``n_iter_`` under its default
    (closed-form) solver, unlike a fitted Lasso model.
    """
    X_train, _X_test, y_train, _y_test = load_train_test()
    ridge = Ridge(alpha=1.0).fit(X_train, y_train)
    lasso = Lasso(alpha=1.0, max_iter=50000).fit(X_train, y_train)
    return {
        "ridge_solver": ridge.solver,
        "ridge_has_n_iter": hasattr(ridge, "n_iter_") and ridge.n_iter_ is not None,
        "lasso_has_n_iter": hasattr(lasso, "n_iter_") and lasso.n_iter_ is not None,
    }


# --------------------------------------------------------------------------
# 8. The corner, in the smallest case that shows it: two correlated features
# --------------------------------------------------------------------------


def two_feature_dataset(n: int = 200, seed: int = 3):
    """Two strongly correlated features and a real, equal-weighted signal,
    small enough to draw. Returns ``(X, y, correlation)``.
    """
    rng = np.random.default_rng(seed)
    x1 = rng.normal(size=n)
    x2 = 0.9 * x1 + rng.normal(scale=0.3, size=n)
    X = np.column_stack([x1, x2])
    y = 2.0 * x1 + 2.0 * x2 + rng.normal(scale=1.0, size=n)
    correlation = float(np.corrcoef(x1, x2)[0, 1])
    return X, y, correlation


def two_feature_corner_demo(alphas):
    """Ridge and lasso coefficients on the two-feature dataset, at each
    alpha, alongside the unregularised (OLS) coefficients for reference.

    Returns ``(ols_coef, rows)`` where each row is
    ``(alpha, ridge_coefs, lasso_coefs)``, coefficients rounded to 4 places.
    """
    from sklearn.linear_model import LinearRegression

    X, y, _correlation = two_feature_dataset()
    ols = LinearRegression().fit(X, y).coef_.round(4).tolist()
    rows = []
    for alpha in alphas:
        ridge = Ridge(alpha=alpha).fit(X, y)
        lasso = Lasso(alpha=alpha, max_iter=50000).fit(X, y)
        rows.append((alpha, ridge.coef_.round(4).tolist(), lasso.coef_.round(4).tolist()))
    return ols, rows

"""Fourteen exercises in what a penalty actually does.

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.
`regularization_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 regularization_lib as r  # noqa: F401  (you will need it)

ALPHA_GRID = [0.001, 0.01, 0.1, 1.0]
PATH_ALPHAS = np.logspace(-3, 2, 60)


def test_01_ridge_never_zeros_lasso_progressively_zeros():
    pytest.skip(
        "Call r.zero_counts_and_r2(ALPHA_GRID) and assert it equals the four "
        "rows in expected-output/measured-values.txt, from "
        "(0.001, 0, 0.3588, 0, 0.3586) to (1.0, 8, 0.2782, 0, 0.357). Then "
        "assert the ridge-zeros column is [0, 0, 0, 0] and the lasso-zeros "
        "column never decreases as alpha grows. Ridge shrinks; lasso "
        "shrinks and selects."
    )


def test_01b_lasso_cv_picks_alpha_and_six_features():
    pytest.skip(
        "Call r.lasso_cv_selection() and assert alpha rounds to 0.07874, "
        "zeros == 4, kept == ['sex', 'bmi', 'bp', 's1', 's3', 's5'], and "
        "r2 == 0.3562. Cross-validation, not eyeballing a curve, is how "
        "alpha actually gets chosen in practice."
    )


def test_02_the_coefficient_path_lasso_hits_exact_zeros_ridge_never_does():
    pytest.skip(
        "Call r.alpha_where_each_lasso_coefficient_first_hits_zero(PATH_ALPHAS), "
        "which returns (zero_at, ridge_ever_zero). Assert every one of the "
        "ten lasso coefficients has a non-None zero_at value, and assert "
        "ridge_ever_zero is False across the same 60-point sweep. A circle "
        "has no corners; a diamond has one on every axis."
    )


def test_02b_the_weakest_coefficients_zero_out_first():
    pytest.skip(
        "From the same zero_at dict, assert the feature with the smallest "
        "zero_at value is 's3' (it rounds to alpha=0.0032) and the feature "
        "with the largest is 'bmi' (alpha=2.4538). The coefficients lasso "
        "can least afford to shrink are the ones it keeps longest."
    )


def test_03_lasso_recovers_the_right_features_at_low_noise():
    pytest.skip(
        "Call r.sparse_recovery(alpha=1.0, noise=1.0, seed=0) against a "
        "known sparse ground truth (5 informative of 20 features) and "
        "assert precision == 1.0, recall == 1.0, n_selected == 5. Repeat at "
        "alpha=0.1, noise=0.1 and assert the same three values. At low "
        "noise, lasso finds exactly the right features -- not merely some "
        "features that predict well."
    )


def test_03b_recovery_degrades_with_noise_and_too_much_penalty():
    pytest.skip(
        "Call r.sparse_recovery(alpha=80.0, noise=30.0, seed=0) and assert "
        "recall == 0.2, n_selected == 1. Call it again at noise=10.0 and "
        "assert recall == 0.0, n_selected == 0 -- too much penalty on noisy "
        "data can zero out the truth entirely. Then call "
        "r.sparse_recovery_across_seeds at (alpha=1.0, noise=1.0) and at "
        "(alpha=1.0, noise=10.0), and assert the second call's mean "
        "precision (0.6792) is lower than the first's (1.0). This is not "
        "one unlucky seed."
    )


def test_04_regularization_requires_scaled_features():
    pytest.skip(
        "Call r.scale_dependence(alpha=1.0). Assert the 'raw' result keeps "
        "all 10 features, 'standardized' keeps 7 "
        "(['sex','bmi','bp','s1','s3','s5','s6']), and 'sklearn_unit_norm' "
        "keeps 3 (['bmi','bp','s5']) -- three different answers, same data, "
        "same alpha. The penalty is applied in whatever units the "
        "coefficients happen to be in, and 'scaled' is not even one "
        "convention: unit-variance and unit-norm disagree too."
    )


def test_05_elasticnet_interpolates_between_ridge_and_lasso():
    pytest.skip(
        "Call r.elasticnet_sweep(alpha=0.1, l1_ratios=[0.0, 0.1, 0.3, 0.5, "
        "0.7, 0.9, 1.0]) and assert it equals the captured rows, ending "
        "(1.0, 3, 0.355). Then assert that row's zero count matches plain "
        "Lasso's own zero count at the same alpha "
        "(r.zero_counts_and_r2([0.1])[0]). l1_ratio=1.0 is not an "
        "approximation to lasso; it is lasso."
    )


def test_05b_ridge_and_elasticnet_do_not_share_an_alpha_scale():
    pytest.skip(
        "Call r.ridge_elasticnet_equivalence(alpha=0.1), which fits "
        "Ridge at alpha * n_train against ElasticNet(alpha=0.1, "
        "l1_ratio=0.0) directly. Assert max_diff < 0.001. Ridge's objective "
        "sums the squared error; ElasticNet's averages it over n_samples -- "
        "so the same alpha means two different penalty strengths until you "
        "correct for it."
    )


def test_06_ridge_splits_the_weight_between_near_duplicates():
    pytest.skip(
        "Get correlation from r.near_duplicate_dataset() and assert it "
        "exceeds 0.999. Call r.ridge_vs_lasso_on_duplicates([1.0]) and, "
        "from the one row, assert the two ridge coefficients differ by "
        "less than 0.15 and sum to within 0.2 of 6.0 -- ridge splits the "
        "true combined weight roughly evenly. Assert lasso's second "
        "coefficient is exactly 0.0 and its first exceeds 5.0 -- lasso "
        "picks one and drops the other."
    )


def test_06b_enough_penalty_zeros_both_duplicates_in_lasso_but_ridge_still_splits():
    pytest.skip(
        "Call r.ridge_vs_lasso_on_duplicates([10.0]). Assert both lasso "
        "coefficients are exactly 0.0. Assert both ridge coefficients "
        "still exceed 2.5 and still differ from each other by less than "
        "0.15. Ridge never produces the corner; more penalty just shrinks "
        "both halves together."
    )


def test_07_ridge_has_a_closed_form_lasso_needs_iterations():
    pytest.skip(
        "Call r.ridge_has_no_iteration_count() and assert "
        "ridge_has_n_iter is False and lasso_has_n_iter is True. Call "
        "r.lasso_iteration_counts(ALPHA_GRID) and assert it equals "
        "{0.001: 368, 0.01: 62, 0.1: 135, 1.0: 6}. Ridge is one "
        "linear-algebra call; lasso is solved iteratively because its "
        "penalty is not differentiable at zero."
    )


def test_08_the_corner_two_correlated_features_and_a_lasso_zero():
    pytest.skip(
        "Call r.two_feature_corner_demo([0.001, 0.5, 1.0, 3.0, 8.0]). At "
        "alpha=3.0, assert lasso's coefficients equal [0.8919, 0.0] exactly "
        "-- it has landed on the axis -- while ridge's second coefficient "
        "still exceeds 1.5 in magnitude. This is the geometry: a diamond's "
        "corner sits on an axis; a circle's does not."
    )


def test_08b_at_a_tiny_alpha_both_models_agree_with_ols():
    pytest.skip(
        "From the same two_feature_corner_demo call (use alphas [0.001, "
        "8.0]), assert both ridge's and lasso's coefficients at alpha=0.001 "
        "are within 0.01 of the OLS coefficients [1.9564, 1.9381]. Then "
        "assert that at alpha=8.0, lasso's coefficients are both exactly "
        "0.0 while ridge's are both still nonzero. As alpha shrinks toward "
        "zero, both penalties vanish and agree with plain least squares."
    )

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

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

import numpy as np

import regularization_lib as r


def test_the_diabetes_split_is_reproducible_and_shaped_right():
    a = r.load_train_test()
    b = r.load_train_test()
    X_train_a, X_test_a, y_train_a, y_test_a = a
    X_train_b, _X_test_b, _y_train_b, _y_test_b = b
    assert X_train_a.shape == (331, 10)
    assert X_test_a.shape == (111, 10)
    assert y_train_a.shape == (331,)
    assert y_test_a.shape == (111,)
    # same seed, same split -- every downstream measurement depends on this
    assert np.array_equal(X_train_a, X_train_b)


def test_the_near_duplicate_columns_really_are_near_duplicate():
    X, y, correlation = r.near_duplicate_dataset()
    assert X.shape == (300, 3)
    assert y.shape == (300,)
    assert correlation > 0.999
    # the third column is not part of the near-duplicate pair
    assert abs(float(np.corrcoef(X[:, 0], X[:, 2])[0, 1])) < 0.3


def test_the_synthetic_sparse_dataset_really_has_five_informative_features():
    from sklearn.datasets import make_regression

    X, y, coef = make_regression(
        n_samples=200, n_features=20, n_informative=5, noise=1.0, coef=True, random_state=0
    )
    assert X.shape == (200, 20)
    assert int(np.sum(coef != 0)) == 5


def test_ridge_solves_directly_and_lasso_iterates():
    info = r.ridge_has_no_iteration_count()
    assert info["ridge_has_n_iter"] is False
    assert info["lasso_has_n_iter"] is True
    counts = r.lasso_iteration_counts([0.001, 1.0])
    # every count is a positive number of iterations, and none of them
    # hit the ceiling -- so no ConvergenceWarning was silently swallowed
    for alpha, n_iter in counts.items():
        assert 0 < n_iter < 50000, f"alpha={alpha} iterated {n_iter} times"

#!/usr/bin/env bash
# Day 151 lab harness: "What the Penalty Does"
#
# 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 manual pytest invocation 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 regularization_lib as r

ALPHA_GRID = [0.001, 0.01, 0.1, 1.0]
PATH_ALPHAS = np.logspace(-3, 2, 60)

errors = []


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


# 1. Ridge never zeros; lasso zeros progressively more
rows = r.zero_counts_and_r2(ALPHA_GRID)
expect(
    "zero counts and r2",
    rows,
    [
        (0.001, 0, 0.3588, 0, 0.3586),
        (0.01, 1, 0.3541, 0, 0.3567),
        (0.1, 3, 0.355, 0, 0.369),
        (1.0, 8, 0.2782, 0, 0.357),
    ],
)
ridge_zeros = [row[3] for row in rows]
lasso_zeros = [row[1] for row in rows]
expect("ridge zeros at every alpha", ridge_zeros, [0, 0, 0, 0])
if not all(a <= b for a, b in zip(lasso_zeros, lasso_zeros[1:])):
    errors.append("lasso zero count decreased somewhere as alpha grew")

# 1b. LassoCV
cv = r.lasso_cv_selection()
expect("LassoCV alpha", round(cv["alpha"], 5), 0.07874)
expect("LassoCV zeros", cv["zeros"], 4)
expect("LassoCV kept", cv["kept"], ["sex", "bmi", "bp", "s1", "s3", "s5"])
expect("LassoCV r2", cv["r2"], 0.3562)

# 2. Coefficient path
zero_at, ridge_ever_zero = r.alpha_where_each_lasso_coefficient_first_hits_zero(PATH_ALPHAS)
if any(v is None for v in zero_at.values()):
    errors.append("not every lasso coefficient zeroed within the sweep")
expect("ridge ever zero", ridge_ever_zero, False)
weakest_first = min(zero_at, key=zero_at.get)
strongest_last = max(zero_at, key=zero_at.get)
expect("weakest coefficient to zero", weakest_first, "s3")
expect("strongest coefficient to zero", strongest_last, "bmi")
expect("s3 zeros at", round(zero_at["s3"], 4), 0.0032)
expect("bmi zeros at", round(zero_at["bmi"], 4), 2.4538)

# 3. Sparse recovery
p1, r1, n1 = r.sparse_recovery(alpha=1.0, noise=1.0, seed=0)
expect("recovery precision, low noise", p1, 1.0)
expect("recovery recall, low noise", r1, 1.0)
expect("recovery n_selected, low noise", n1, 5)

# 3b. Recovery degrades
_p, r_high, n_high = r.sparse_recovery(alpha=80.0, noise=30.0, seed=0)
expect("recall, heavy penalty + noise", r_high, 0.2)
expect("n_selected, heavy penalty + noise", n_high, 1)
_p2, r_none, n_none = r.sparse_recovery(alpha=80.0, noise=10.0, seed=0)
expect("recall, heavy penalty zeroes truth", r_none, 0.0)
expect("n_selected, heavy penalty zeroes truth", n_none, 0)
mp_low, mr_low = r.sparse_recovery_across_seeds(alpha=1.0, noise=1.0)
mp_high, mr_high = r.sparse_recovery_across_seeds(alpha=1.0, noise=10.0)
expect("mean precision, low noise, 10 seeds", mp_low, 1.0)
expect("mean precision, high noise, 10 seeds", mp_high, 0.6792)
if not (mp_high < mp_low):
    errors.append("mean precision did not fall with more noise")

# 4. Scale dependence
scale = r.scale_dependence(alpha=1.0)
expect("raw n_kept", scale["raw"]["n_kept"], 10)
expect("standardized n_kept", scale["standardized"]["n_kept"], 7)
expect("sklearn_unit_norm n_kept", scale["sklearn_unit_norm"]["n_kept"], 3)
expect(
    "standardized kept set",
    scale["standardized"]["kept"],
    ["sex", "bmi", "bp", "s1", "s3", "s5", "s6"],
)
expect("sklearn_unit_norm kept set", scale["sklearn_unit_norm"]["kept"], ["bmi", "bp", "s5"])

# 5. ElasticNet
en_rows = r.elasticnet_sweep(alpha=0.1, l1_ratios=[0.0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.0])
expect(
    "elasticnet sweep",
    en_rows,
    [
        (0.0, 0, 0.0555),
        (0.1, 0, 0.0605),
        (0.3, 0, 0.0741),
        (0.5, 0, 0.0963),
        (0.7, 1, 0.1389),
        (0.9, 0, 0.2511),
        (1.0, 3, 0.355),
    ],
)

# 5b. Ridge/ElasticNet alpha scale
ridge_head, elastic_head, max_diff = r.ridge_elasticnet_equivalence(alpha=0.1)
if not (max_diff < 0.001):
    errors.append(f"ridge/elasticnet correction did not converge: max_diff={max_diff}")
expect("ridge head", ridge_head, [6.3098, 1.229, 22.6076])
expect("elastic head", elastic_head, [6.3097, 1.229, 22.6076])

# 6. Near duplicates
_X, _y, correlation = r.near_duplicate_dataset()
if not (correlation > 0.999):
    errors.append(f"near-duplicate correlation too low: {correlation}")
dup_rows_1 = r.ridge_vs_lasso_on_duplicates([1.0])
_alpha, ridge_c, lasso_c = dup_rows_1[0]
if abs(ridge_c[0] - ridge_c[1]) >= 0.15:
    errors.append("ridge did not split the duplicate weight evenly at alpha=1.0")
if lasso_c[1] != 0.0:
    errors.append("lasso did not zero the second duplicate at alpha=1.0")

# 6b. High alpha
dup_rows_10 = r.ridge_vs_lasso_on_duplicates([10.0])
_alpha10, ridge_c10, lasso_c10 = dup_rows_10[0]
expect("lasso zeros both duplicates at alpha=10", lasso_c10, [0.0, 0.0, 0.0])
if not (ridge_c10[0] > 2.5 and ridge_c10[1] > 2.5):
    errors.append("ridge zeroed a duplicate at alpha=10, which it must never do")

# 7. Closed form vs iterative
info = r.ridge_has_no_iteration_count()
expect("ridge has n_iter_", info["ridge_has_n_iter"], False)
expect("lasso has n_iter_", info["lasso_has_n_iter"], True)
counts = r.lasso_iteration_counts(ALPHA_GRID)
expect("lasso iteration counts", counts, {0.001: 368, 0.01: 62, 0.1: 135, 1.0: 6})

# 8. The corner
ols, corner_rows = r.two_feature_corner_demo([0.001, 0.5, 1.0, 3.0, 8.0])
expect("OLS reference coefficients", ols, [1.9564, 1.9381])
by_alpha = {a: (rc, lc) for a, rc, lc in corner_rows}
ridge3, lasso3 = by_alpha[3.0]
expect("lasso lands on the axis at alpha=3.0", lasso3, [0.8919, 0.0])
if ridge3[0] == 0.0 or ridge3[1] == 0.0:
    errors.append("ridge produced an exact zero, which it must never do")

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 regularization_lib, no pytest involved"
else
  fail "direct library checks failed"
  echo "$DIRECT_CHECK" | sed 's/^/    /'
fi

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

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

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

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

echo ""
echo "7. Proof the harness can fail"
SCRATCH=$(mktemp -d "${TMPDIR:-/tmp}/d151-scratch.XXXXXX")
cp examples/*.py "$SCRATCH"/
SCRATCH_OUT=$("$PYTEST" "$SCRATCH" -q 2>&1)
if echo "$SCRATCH_OUT" | tail -1 | grep -qE "^18 passed"; then
  ok "scratch copy of examples/ passes before it is broken"
else
  fail "scratch copy did not pass before being broken: $(echo "$SCRATCH_OUT" | tail -3)"
fi
"$PYTHON" - "$SCRATCH/test_regularization_claims.py" <<'PYEOF'
import sys
path = sys.argv[1]
text = open(path).read()
needle = "assert lasso_x2 == 0.0"
replacement = "assert lasso_x2 == 999.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_06_ridge_splits_the_weight_between_near_duplicates"; then
  ok "breaking exercise 6's assertion produces a non-zero exit and names the failing test"
else
  fail "broken copy did not fail as expected (exit=$BROKEN_STATUS)"
fi
rm -rf "$SCRATCH"

echo ""
echo "8. The direction of every result holds beyond the quoted alpha or seed"
DIRECTION=$("$PYTHON" - <<'PYEOF'
import sys
sys.path.insert(0, "examples")
import numpy as np
import regularization_lib as r

problems = []

# Ridge never zeros a coefficient, at alphas the lesson does not quote.
_lasso_path, ridge_path = r.coefficient_path(np.logspace(-4, 3, 25))
if np.any(ridge_path == 0.0):
    problems.append("ridge produced an exact zero at an unquoted alpha")

# Lasso's zero count is monotone non-decreasing over a different, denser grid.
alphas = np.logspace(-3, 0.5, 12)
lasso_zeros = []
X_train, X_test, y_train, y_test = r.load_train_test()
from sklearn.linear_model import Lasso
for a in alphas:
    lasso_zeros.append(int(np.sum(Lasso(alpha=a, max_iter=50000).fit(X_train, y_train).coef_ == 0)))
if not all(a <= b for a, b in zip(lasso_zeros, lasso_zeros[1:])):
    problems.append(f"lasso zero count was not monotone over a denser alpha grid: {lasso_zeros}")

# Near-duplicate splitting holds at other dataset seeds too.
for seed in (1, 2, 3):
    X, y, corr = r.near_duplicate_dataset(seed=seed)
    if corr <= 0.999:
        problems.append(f"seed {seed}: near-duplicate correlation too low ({corr})")
    from sklearn.linear_model import Ridge
    ridge = Ridge(alpha=1.0).fit(X, y)
    lasso = Lasso(alpha=1.0, max_iter=50000).fit(X, y)
    if abs(ridge.coef_[0] - ridge.coef_[1]) >= 0.2:
        problems.append(f"seed {seed}: ridge did not split the duplicate weight")
    # lasso need not land on an EXACT zero at every seed -- these columns
    # are correlated at 0.9999, not identically collinear -- but it must
    # produce a heavily asymmetric split rather than ridge's near-even one
    small, large = sorted(abs(c) for c in lasso.coef_[:2])
    if large < 20 * max(small, 1e-6):
        problems.append(f"seed {seed}: lasso did not produce an asymmetric split ({lasso.coef_[:2]})")

if problems:
    for p in problems:
        print("ERROR:", p)
else:
    print("every direction held")
PYEOF
)
if [ "$DIRECTION" = "every direction held" ]; then
  ok "ridge never-zeros, lasso monotone zeroing, and duplicate splitting hold beyond the quoted alphas and seeds"
else
  fail "a direction failed beyond the quoted alpha or 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