Hands-on lab — Day 153: Linear Regression from Scratch

"""Ordinary least squares, three ways: what the library was doing for you.

Three implementations of the same fit -- the normal equations, an
lstsq-based solve, and gradient descent -- plus a scikit-learn-compatible
estimator wrapping all three, and then a direct measurement of where they
agree and where they do not.

Day 111 already derived gradient descent, its update rule, the stability
condition ``|1 - eta * a| < 1`` for an eigenvalue ``a`` of the loss's
Hessian, and the condition number as the ratio of the Hessian's largest and
smallest eigenvalues. None of that is re-derived here; it is applied.

Day 146 already established the scikit-learn estimator API contract --
``fit``/``predict``, learned attributes with a trailing underscore -- and
measured that a from-scratch estimator raises ``AttributeError`` on
``__sklearn_tags__`` inside ``Pipeline.predict()`` and ``cross_val_score``
unless it inherits ``BaseEstimator``. ``OLSRegressor`` below inherits it for
exactly that reason, without re-teaching the contract.

Everything here is deterministic given a seed.
"""

from __future__ import annotations

import numpy as np

from sklearn.base import BaseEstimator, RegressorMixin
from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression
from sklearn.utils.validation import check_array, check_is_fitted, check_X_y

# --------------------------------------------------------------------------
# 1. Loading the data, and the two condition numbers
# --------------------------------------------------------------------------


def load_diabetes_data(scaled: bool = True):
    """The bundled diabetes regression dataset: X (442, 10), y (442,).

    ``scaled=True`` (the default) gives sklearn's own mean-centred,
    unit-L2-norm-scaled columns. ``scaled=False`` gives the raw clinical
    units -- age in years, sex coded 1/2, six serum measurements on very
    different scales -- which is the badly scaled version used to show
    gradient descent struggling.
    """
    data = load_diabetes(scaled=scaled)
    return data.data, data.target


def add_intercept_column(X: np.ndarray) -> np.ndarray:
    """Prepend a column of ones, the classic way to fold in an intercept."""
    ones = np.ones((X.shape[0], 1))
    return np.hstack([ones, X])


def condition_numbers(A: np.ndarray) -> tuple[float, float]:
    """``(cond(A), cond(A'A))`` -- the second is exactly the square of the first.

    ``cond`` here is the ratio of largest to smallest singular value. Squaring
    a matrix's condition number is the textbook reason the normal equations
    lose precision that a direct solve of ``A`` does not.
    """
    cond_a = float(np.linalg.cond(A))
    cond_ata = float(np.linalg.cond(A.T @ A))
    return cond_a, cond_ata


# --------------------------------------------------------------------------
# 2. Three ways to fit: normal equations, lstsq, and sklearn as the referee
# --------------------------------------------------------------------------


def fit_normal_equations(A: np.ndarray, y: np.ndarray) -> np.ndarray:
    """Solve the normal equations ``A'A b = A'y`` directly.

    This is the formula from the textbook page: form ``A'A``, form ``A'y``,
    solve the square system. It is also, as this module measures, the least
    numerically careful of the three routes here.
    """
    return np.linalg.solve(A.T @ A, A.T @ y)


def fit_lstsq(A: np.ndarray, y: np.ndarray) -> np.ndarray:
    """Solve the same least-squares problem via ``numpy.linalg.lstsq``.

    ``lstsq`` factors ``A`` directly (an SVD-based solve under the hood)
    rather than forming and inverting ``A'A``, which is why it does not
    inherit the squared condition number of the normal equations.
    """
    coef, *_ = np.linalg.lstsq(A, y, rcond=None)
    return coef


def sklearn_reference_fit(X: np.ndarray, y: np.ndarray) -> np.ndarray:
    """Fit scikit-learn's own ``LinearRegression`` and return ``[intercept, *coef]``.

    Used as the referee throughout this module: not because it is a fourth
    algorithm, but because it is the thing every from-scratch attempt here is
    being measured against.
    """
    model = LinearRegression().fit(X, y)
    return np.concatenate([[float(model.intercept_)], model.coef_])


def max_abs_difference(a: np.ndarray, b: np.ndarray) -> float:
    return float(np.max(np.abs(a - b)))


# --------------------------------------------------------------------------
# 3. The dramatic case: a near-duplicate column
# --------------------------------------------------------------------------


def make_dramatic_collinear_dataset(n: int = 100, seed: int = 0, noise_scale: float = 1e-7):
    """Three random predictors plus a fourth that is almost the first.

    The fourth column is column 0 plus a sliver of noise -- 1e-7 in scale,
    far smaller than any real measurement error. The true coefficients are
    [1, 2, 3, 4]; a well-conditioned fit should recover something close to
    that. It does not, for two of the three methods below.
    """
    rng = np.random.default_rng(seed)
    X = rng.normal(size=(n, 3))
    duplicate = X[:, 0] + rng.normal(scale=noise_scale, size=n)
    X = np.column_stack([X, duplicate])
    true_coef = np.array([1.0, 2.0, 3.0, 4.0])
    y = X @ true_coef + rng.normal(scale=0.1, size=n)
    return X, y, true_coef


# --------------------------------------------------------------------------
# 4. Gradient descent, and the Day 111 stability threshold arriving here
# --------------------------------------------------------------------------


def standardize(X: np.ndarray) -> np.ndarray:
    """Zero mean, unit standard deviation, column by column."""
    mu = X.mean(axis=0)
    sd = X.std(axis=0)
    return (X - mu) / sd


def center(X: np.ndarray, y: np.ndarray):
    """Subtract each column's mean from X and y's mean from y.

    Centring removes the need for a separate intercept term: the fitted
    hyperplane through centred data passes through the origin, and the
    intercept is recovered afterwards as ``y.mean() - X.mean(axis=0) @ coef``.
    """
    X_mean = X.mean(axis=0)
    y_mean = y.mean()
    return X - X_mean, y - y_mean, X_mean, y_mean


def hessian_eigenvalues(X: np.ndarray, n: int) -> np.ndarray:
    """Eigenvalues of the mean-squared-error loss's Hessian, ``(2/n) X'X``.

    Day 111 established that gradient descent on a quadratic loss is stable
    exactly when ``|1 - eta * a| < 1`` for every eigenvalue ``a`` of this
    Hessian, and that the ratio of its largest to smallest eigenvalue is the
    condition number governing how slowly the slowest direction converges.
    """
    XtX = X.T @ X
    return np.linalg.eigvalsh((2.0 / n) * XtX)


def stability_threshold(X: np.ndarray) -> float:
    """The largest learning rate for which every eigenvalue keeps |1 - eta*a| < 1.

    From Day 111's condition: stability requires ``eta < 2 / a`` for every
    eigenvalue ``a``, so the binding constraint is the largest eigenvalue.
    """
    n = X.shape[0]
    eig_max = float(hessian_eigenvalues(X, n).max())
    return 2.0 / eig_max


def fit_gradient_descent(X: np.ndarray, y: np.ndarray, lr: float, n_iter: int) -> np.ndarray:
    """Plain batch gradient descent on the centred mean-squared-error loss.

    Starts at the zero vector. ``X`` and ``y`` are assumed already centred,
    so no intercept term appears in this loop -- see ``center`` above.
    """
    n, p = X.shape
    coef = np.zeros(p)
    for _ in range(n_iter):
        grad = (2.0 / n) * (X.T @ (X @ coef - y))
        coef = coef - lr * grad
    return coef


def iters_to_tolerance(X: np.ndarray, y: np.ndarray, lr: float, target: np.ndarray, tol: float, max_iter: int):
    """How many gradient-descent iterations until every coefficient is within tol of target.

    Returns ``(iterations, final_coef)``, or ``(None, final_coef)`` if the
    tolerance was never reached within ``max_iter`` -- either because
    convergence needs more steps, or because the run diverged. A run that
    stops producing finite numbers returns ``("diverged", coef)``.
    """
    n, p = X.shape
    coef = np.zeros(p)
    for i in range(1, max_iter + 1):
        grad = (2.0 / n) * (X.T @ (X @ coef - y))
        coef = coef - lr * grad
        if not np.all(np.isfinite(coef)):
            return "diverged", coef
        if float(np.max(np.abs(coef - target))) < tol:
            return i, coef
    return None, coef


# --------------------------------------------------------------------------
# 5. Counting operations instead of timing them
# --------------------------------------------------------------------------


def normal_equation_op_count(n: int, p: int) -> int:
    """Multiply-adds to form A'A (n*p^2) plus solve the p-by-p system (p^3)."""
    return n * p * p + p**3


def gradient_descent_op_count(n: int, p: int, iterations: int) -> int:
    """Multiply-adds for `iterations` steps: two n*p matrix-vector products each."""
    return 2 * n * p * iterations


# --------------------------------------------------------------------------
# 6. A scikit-learn-compatible estimator
# --------------------------------------------------------------------------


class OLSRegressor(RegressorMixin, BaseEstimator):
    """Ordinary least squares, fit by one of three methods, as a real estimator.

    Inherits ``RegressorMixin`` and ``BaseEstimator`` because Day 146 already
    measured what happens without them: ``fit``, ``predict`` and ``score``
    work perfectly well called directly, but ``Pipeline.predict()`` and
    ``cross_val_score`` both raise ``AttributeError`` on ``__sklearn_tags__``,
    which only ``BaseEstimator`` supplies. That lesson is assumed, not
    repeated.

    ``fit_intercept=True`` centres X and y, fits the centred problem, and
    recovers the intercept afterwards -- rather than appending a column of
    ones -- because the two are measured elsewhere in this module to agree
    to about ten decimal places and centring avoids growing the design
    matrix by one column.
    """

    def __init__(self, method: str = "lstsq", fit_intercept: bool = True, lr: float = 0.1, n_iter: int = 1000):
        self.method = method
        self.fit_intercept = fit_intercept
        self.lr = lr
        self.n_iter = n_iter

    def fit(self, X, y):
        X, y = check_X_y(X, y)
        if self.fit_intercept:
            Xc, yc, X_mean, y_mean = center(X, y)
        else:
            Xc, yc = X, y
            X_mean, y_mean = np.zeros(X.shape[1]), 0.0

        if self.method == "normal":
            coef = fit_normal_equations(Xc, yc)
        elif self.method == "lstsq":
            coef = fit_lstsq(Xc, yc)
        elif self.method == "gd":
            coef = fit_gradient_descent(Xc, yc, self.lr, self.n_iter)
        else:
            raise ValueError(f"unknown method {self.method!r}; use 'normal', 'lstsq' or 'gd'")

        self.coef_ = coef
        self.intercept_ = float(y_mean - X_mean @ coef) if self.fit_intercept else 0.0
        self.n_features_in_ = X.shape[1]
        return self

    def predict(self, X):
        check_is_fitted(self, "coef_")
        X = check_array(X)
        return X @ self.coef_ + self.intercept_


def run_check_estimator(estimator):
    """Run scikit-learn's estimator_checks suite and tally results by name.

    Returns ``(passed, failed, skipped)`` where ``failed`` and ``skipped``
    are lists of ``(check_name, message)`` pairs -- never silently discarded.
    """
    from sklearn.utils.estimator_checks import check_estimator

    results = []

    def record(*, estimator, check_name, exception, status, expected_to_fail, expected_to_fail_reason):
        results.append((check_name, status, str(exception)[:200] if exception else None))

    check_estimator(estimator, on_fail=None, on_skip=None, callback=record)
    passed = [name for name, status, _msg in results if status == "passed"]
    failed = [(name, msg) for name, status, msg in results if status == "failed"]
    skipped = [(name, msg) for name, status, msg in results if status == "skipped"]
    return passed, failed, skipped


# --------------------------------------------------------------------------
# 7. fit_intercept two ways: centring versus an appended column
# --------------------------------------------------------------------------


def fit_intercept_two_ways(X: np.ndarray, y: np.ndarray):
    """Compare centring against appending a ones column, on the same data.

    Returns ``(coef_column, intercept_column, coef_centred, intercept_centred)``.
    """
    n = X.shape[0]
    A = add_intercept_column(X)
    beta_column = fit_normal_equations(A, y)
    intercept_column, coef_column = float(beta_column[0]), beta_column[1:]

    Xc, yc, X_mean, y_mean = center(X, y)
    coef_centred = fit_normal_equations(Xc, yc)
    intercept_centred = float(y_mean - X_mean @ coef_centred)

    return coef_column, intercept_column, coef_centred, intercept_centred

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

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

import sys
from pathlib import Path

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

import numpy as np  # noqa: E402
from sklearn.linear_model import LinearRegression  # 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 153 -- linear regression from scratch, measured")
    print("=" * 52)

    X, y = r.load_diabetes_data(scaled=True)
    n, p = X.shape

    rule("1. Three closed forms on well-conditioned data")
    A = r.add_intercept_column(X)
    beta_ne = r.fit_normal_equations(A, y)
    beta_lstsq = r.fit_lstsq(A, y)
    beta_sk = r.sklearn_reference_fit(X, y)
    gap_ne = r.max_abs_difference(beta_ne, beta_sk)
    gap_lstsq = r.max_abs_difference(beta_lstsq, beta_sk)
    print(f"  max |normal equations - sklearn| : {gap_ne:.4e}")
    print(f"  max |lstsq            - sklearn| : {gap_lstsq:.4e}")
    print(f"  lstsq is {gap_ne / gap_lstsq:.1f}x closer to sklearn's own answer")

    rule("1b. The normal equations square the condition number")
    cond_a, cond_ata = r.condition_numbers(A)
    print(f"  condition number of X (with intercept) : {cond_a:.4f}")
    print(f"  condition number of X'X                : {cond_ata:.4f}")
    print(f"  cond(X'X) / cond(X)^2                   : {cond_ata / cond_a**2:.10f}")

    rule("2. A near-duplicate column, and three very different answers")
    Xd, yd, true_coef = r.make_dramatic_collinear_dataset(n=100, seed=0)
    Ad = r.add_intercept_column(Xd)
    beta_ne_d = r.fit_normal_equations(Ad, yd)
    beta_lstsq_d = r.fit_lstsq(Ad, yd)
    beta_sk_d = r.sklearn_reference_fit(Xd, yd)
    print("  true coefficients (intercept, c0, c1, c2, c3-duplicate)")
    print(f"    = [0.0, {true_coef[0]:.1f}, {true_coef[1]:.1f}, {true_coef[2]:.1f}, {true_coef[3]:.1f}]")
    print(f"  normal-equation coefficients : {np.round(beta_ne_d, 3).tolist()}")
    print(f"  lstsq coefficients           : {np.round(beta_lstsq_d, 3).tolist()}")
    print(f"  sklearn coefficients         : {np.round(beta_sk_d, 3).tolist()}")

    rule("2b. And the squaring relationship itself becomes hard to verify")
    cond_ad, cond_atad = r.condition_numbers(Ad)
    print(f"  condition number of X (with intercept) : {cond_ad:.4e}")
    print(f"  condition number of X'X                : {cond_atad:.4e}")
    print(f"  cond(X'X) / cond(X)^2                   : {cond_atad / cond_ad**2:.4f}")

    rule("3. Gradient descent versus the closed form, standardized features")
    Xs = r.standardize(X)
    yc = y - y.mean()
    target = r.fit_normal_equations(Xs, yc)
    threshold = r.stability_threshold(Xs)
    eig_s = r.hessian_eigenvalues(Xs, n)
    print(f"  stability threshold (Day 111, |1 - eta*a| < 1) : {threshold:.4f}")
    print(f"  Hessian eigenvalue ratio (max/min)             : {eig_s.max() / eig_s.min():.4f}")
    lr = 0.2
    iters_3, _ = r.iters_to_tolerance(Xs, yc, lr, target, 5e-4, 200_000)
    iters_6, _ = r.iters_to_tolerance(Xs, yc, lr, target, 5e-7, 200_000)
    iters_9, _ = r.iters_to_tolerance(Xs, yc, lr, target, 5e-10, 200_000)
    print(f"  at lr={lr} (80 percent of threshold):")
    print(f"    iterations to agree to 3 decimals : {iters_3}")
    print(f"    iterations to agree to 6 decimals : {iters_6}")
    print(f"    iterations to agree to 9 decimals : {iters_9}")

    rule("3b. The same setup on raw, unscaled features")
    Xraw, yraw = r.load_diabetes_data(scaled=False)
    Xrc, yrc, _, _ = r.center(Xraw, yraw)
    eig_r = r.hessian_eigenvalues(Xrc, n)
    raw_threshold = r.stability_threshold(Xrc)
    target_r = r.fit_normal_equations(Xrc, yrc)
    print(f"  stability threshold, raw features : {raw_threshold:.6e}")
    print(f"  Hessian eigenvalue ratio, raw      : {eig_r.max() / eig_r.min():.4f}")
    status, coef_r = r.iters_to_tolerance(Xrc, yrc, raw_threshold * 0.95, target_r, 5e-4, 200_000)
    print(f"  at 95 percent of ITS OWN threshold, after 200000 iterations:")
    print(f"    converged to 3 decimals? {status is not None}")
    print(f"    remaining max |coef - closed form| : {r.max_abs_difference(coef_r, target_r):.4f}")

    rule("4. The stability threshold predicts divergence exactly")
    below_status, _ = r.iters_to_tolerance(Xs, yc, threshold * 0.8, target, 1e-9, 20_000)
    above_status, above_coef = r.iters_to_tolerance(Xs, yc, threshold * 1.02, target, 1e-9, 20_000)
    print(f"  at 80 percent of threshold  (lr={threshold * 0.8:.4f}): converges in {below_status} iterations")
    print(f"  at 102 percent of threshold (lr={threshold * 1.02:.4f}): {above_status}, finite={np.all(np.isfinite(above_coef))}")

    rule("5. Operations, not time")
    ops_normal = r.normal_equation_op_count(n, p + 1)
    ops_gd = r.gradient_descent_op_count(n, p, iters_9)
    print(f"  normal-equation operations (form X'X, O(n p^2), plus solve, O(p^3)) : {ops_normal:,}")
    print(f"  gradient-descent operations ({iters_9} iterations, O(n p) each)      : {ops_gd:,}")
    print(f"  ratio                                                              : {ops_gd / ops_normal:.2f}x")

    rule("6. A scikit-learn-compatible estimator, checked by the library itself")
    sk = LinearRegression().fit(Xs, y)
    normal_est = r.OLSRegressor(method="normal").fit(Xs, y)
    lstsq_est = r.OLSRegressor(method="lstsq").fit(Xs, y)
    gd_est = r.OLSRegressor(method="gd", lr=0.2, n_iter=8000).fit(Xs, y)
    print(f"  max |normal-method coef  - sklearn| : {r.max_abs_difference(normal_est.coef_, sk.coef_):.4e}")
    print(f"  max |lstsq-method coef   - sklearn| : {r.max_abs_difference(lstsq_est.coef_, sk.coef_):.4e}")
    print(f"  max |gd-method coef      - sklearn| : {r.max_abs_difference(gd_est.coef_, sk.coef_):.4e}")
    passed, failed, skipped = r.run_check_estimator(r.OLSRegressor())
    print(f"  check_estimator: {len(passed)} passed, {len(failed)} failed, {len(skipped)} skipped, {len(passed) + len(failed) + len(skipped)} total")
    for name, _msg in failed:
        print(f"    FAILED : {name}")
    for name, _msg in skipped:
        print(f"    SKIPPED: {name}")

    rule("7. fit_intercept: centring versus an appended column of ones")
    coef_col, intercept_col, coef_centred, intercept_centred = r.fit_intercept_two_ways(X, y)
    print(f"  max |coef_column - coef_centred| : {r.max_abs_difference(coef_col, coef_centred):.4e}")
    print(f"  |intercept_column - intercept_centred| : {abs(intercept_col - intercept_centred):.4e}")


if __name__ == "__main__":
    main()

"""Ten exercises in what a linear-regression library was doing for you.

Reference solutions. Read `starter/00_brief.md` and
`starter/test_regression_claims.py` for the exercise version, where each of
these bodies is a `pytest.skip` naming exactly what to build.

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 pytest

import regression_lib as r


@pytest.fixture(scope="module")
def diabetes():
    return r.load_diabetes_data(scaled=True)


@pytest.fixture(scope="module")
def diabetes_raw():
    return r.load_diabetes_data(scaled=False)


def test_01_lstsq_agrees_with_sklearn_a_hundred_times_more_closely_than_the_normal_equations(diabetes):
    X, y = diabetes
    A = r.add_intercept_column(X)
    beta_ne = r.fit_normal_equations(A, y)
    beta_lstsq = r.fit_lstsq(A, y)
    beta_sk = r.sklearn_reference_fit(X, y)

    gap_ne = r.max_abs_difference(beta_ne, beta_sk)
    gap_lstsq = r.max_abs_difference(beta_lstsq, beta_sk)

    assert gap_ne == pytest.approx(1.2153e-10, rel=0.05)
    assert gap_lstsq == pytest.approx(1.199e-12, rel=0.05)
    # lstsq is roughly a hundred times closer to sklearn's own answer
    assert 50 < gap_ne / gap_lstsq < 200


def test_01b_the_normal_equations_condition_number_is_exactly_the_square(diabetes):
    X, y = diabetes
    A = r.add_intercept_column(X)
    cond_a, cond_ata = r.condition_numbers(A)

    assert cond_a == pytest.approx(227.2248, rel=1e-4)
    assert cond_ata == pytest.approx(51631.1119, rel=1e-4)
    # cond(A'A) is the square of cond(A), to numerical precision here
    assert cond_ata / cond_a**2 == pytest.approx(1.0, abs=1e-8)


def test_02_a_near_duplicate_column_makes_the_normal_equations_explode():
    X, y, true_coef = r.make_dramatic_collinear_dataset(n=100, seed=0)
    A = r.add_intercept_column(X)
    beta_ne = r.fit_normal_equations(A, y)
    beta_lstsq = r.fit_lstsq(A, y)
    beta_sk = r.sklearn_reference_fit(X, y)

    # the true coefficients for the two duplicated columns are 1 and 4;
    # both from-scratch closed forms split that weight into something
    # unrecognisable, in opposite directions
    assert beta_ne[1] > 1e5
    assert beta_ne[4] < -1e5
    assert beta_lstsq[1] > 1e5
    assert beta_lstsq[4] < -1e5

    # sklearn's own LinearRegression -- an SVD-based minimum-norm solve --
    # stays sane and splits the weight evenly between the near-duplicates
    assert beta_sk[1] == pytest.approx(2.5, abs=0.05)
    assert beta_sk[4] == pytest.approx(2.5, abs=0.05)
    # and it recovers the two non-duplicated coefficients accurately
    assert beta_sk[2] == pytest.approx(2.0, abs=0.05)
    assert beta_sk[3] == pytest.approx(3.0, abs=0.05)


def test_02b_even_the_squaring_relationship_becomes_hard_to_verify_here():
    X, y, _true_coef = r.make_dramatic_collinear_dataset(n=100, seed=0)
    A = r.add_intercept_column(X)
    cond_a, cond_ata = r.condition_numbers(A)

    # several orders of magnitude worse than the diabetes case above
    assert cond_a > 1e6
    assert cond_ata > 1e13

    # the theoretical relationship is exact, but computing the smallest
    # singular value of an already near-singular matrix is itself
    # imprecise, so the measured ratio drifts noticeably from 1.0 -- unlike
    # the clean diabetes case, where it held to twelve decimal places
    ratio = cond_ata / cond_a**2
    assert 0.9 < ratio < 1.0
    assert abs(ratio - 1.0) > 1e-4


def test_03_gradient_descent_reaches_the_closed_form_at_three_six_and_nine_decimals(diabetes):
    X, y = diabetes
    Xs = r.standardize(X)
    yc = y - y.mean()
    target = r.fit_normal_equations(Xs, yc)

    threshold = r.stability_threshold(Xs)
    assert threshold == pytest.approx(0.2485, rel=1e-3)

    learning_rate = 0.2  # about 80 percent of the stability threshold
    iters_3, _ = r.iters_to_tolerance(Xs, yc, learning_rate, target, 5e-4, 200_000)
    iters_6, _ = r.iters_to_tolerance(Xs, yc, learning_rate, target, 5e-7, 200_000)
    iters_9, _ = r.iters_to_tolerance(Xs, yc, learning_rate, target, 5e-10, 200_000)

    assert iters_3 == 3263
    assert iters_6 == 5277
    assert iters_9 == 7291
    # each extra three decimal places costs a comparable number of further
    # iterations, not an exponentially larger one -- linear convergence
    assert iters_6 - iters_3 < iters_3
    assert iters_9 - iters_6 < iters_6


def test_03b_the_same_setup_on_unscaled_features_barely_moves(diabetes, diabetes_raw):
    X, y = diabetes
    Xraw, yraw = diabetes_raw
    Xs = r.standardize(X)

    eig_scaled = r.hessian_eigenvalues(Xs, Xs.shape[0])
    Xrc, yrc, _, _ = r.center(Xraw, yraw)
    eig_raw = r.hessian_eigenvalues(Xrc, Xrc.shape[0])

    ratio_scaled = float(eig_scaled.max() / eig_scaled.min())
    ratio_raw = float(eig_raw.max() / eig_raw.min())

    assert ratio_scaled == pytest.approx(470.08, rel=0.01)
    assert ratio_raw == pytest.approx(76278.96, rel=0.01)
    # the raw features are over a hundred times worse conditioned, in the
    # Hessian-eigenvalue sense Day 111 already established
    assert ratio_raw / ratio_scaled > 100

    raw_threshold = r.stability_threshold(Xrc)
    target_raw = r.fit_normal_equations(Xrc, yrc)
    status, coef = r.iters_to_tolerance(Xrc, yrc, raw_threshold * 0.95, target_raw, 5e-4, 200_000)

    # a learning rate at 95 percent of ITS OWN stability threshold --
    # stable in principle -- still has not reached even one decimal place
    # of agreement after 200,000 iterations, because the slowest direction
    # is governed by the smallest eigenvalue, not the largest
    assert status is None
    assert r.max_abs_difference(coef, target_raw) > 0.1


def test_04_the_day_111_stability_threshold_predicts_divergence_exactly(diabetes):
    X, y = diabetes
    Xs = r.standardize(X)
    yc = y - y.mean()
    target = r.fit_normal_equations(Xs, yc)
    threshold = r.stability_threshold(Xs)

    # comfortably below threshold: converges
    below_status, below_coef = r.iters_to_tolerance(Xs, yc, threshold * 0.8, target, 1e-9, 20_000)
    assert below_status == 7132

    # just past threshold: diverges to non-finite values
    above_status, above_coef = r.iters_to_tolerance(Xs, yc, threshold * 1.02, target, 1e-9, 20_000)
    assert above_status == "diverged"
    assert not np.all(np.isfinite(above_coef))


def test_05_the_closed_form_needs_a_thousand_times_fewer_operations(diabetes):
    X, _y = diabetes
    n, p = X.shape

    ops_normal = r.normal_equation_op_count(n, p + 1)  # +1 for the intercept column
    ops_gd = r.gradient_descent_op_count(n, p, iterations=7291)  # 9-decimal convergence, exercise 3

    assert ops_normal == 54813
    assert ops_gd == 64_452_440
    ratio = ops_gd / ops_normal
    assert 1000 < ratio < 1300


def test_06_the_estimator_matches_sklearn_and_check_estimator_names_two_failures(diabetes):
    from sklearn.linear_model import LinearRegression

    X, y = diabetes
    Xs = r.standardize(X)
    sk = LinearRegression().fit(Xs, y)

    normal_est = r.OLSRegressor(method="normal").fit(Xs, y)
    lstsq_est = r.OLSRegressor(method="lstsq").fit(Xs, y)
    gd_est = r.OLSRegressor(method="gd", lr=0.2, n_iter=8000).fit(Xs, y)

    assert r.max_abs_difference(normal_est.coef_, sk.coef_) < 1e-9
    assert r.max_abs_difference(lstsq_est.coef_, sk.coef_) < 1e-9
    assert r.max_abs_difference(gd_est.coef_, sk.coef_) < 1e-8

    passed, failed, skipped = r.run_check_estimator(r.OLSRegressor())
    failed_names = {name for name, _msg in failed}
    skipped_names = {name for name, _msg in skipped}

    assert len(passed) == 48
    assert failed_names == {"check_n_features_in_after_fitting", "check_dtype_object"}
    assert skipped_names == {"check_array_api_input", "check_regressor_data_not_an_array"}
    assert len(passed) + len(failed) + len(skipped) == 52


def test_07_centring_and_appending_a_column_agree_to_ten_decimal_places(diabetes):
    X, y = diabetes
    coef_col, intercept_col, coef_centred, intercept_centred = r.fit_intercept_two_ways(X, y)

    assert r.max_abs_difference(coef_col, coef_centred) < 1e-9
    assert abs(intercept_col - intercept_centred) < 1e-9

"""Machinery checks -- already solved, in both starter/ and examples/.

These confirm the library behaves as documented. They are not the
exercises; `test_regression_claims.py` is.
"""

import numpy as np

import regression_lib as r


def test_add_intercept_column_prepends_a_column_of_ones():
    X = np.arange(6.0).reshape(3, 2)
    A = r.add_intercept_column(X)
    assert A.shape == (3, 3)
    assert np.array_equal(A[:, 0], np.ones(3))
    assert np.array_equal(A[:, 1:], X)


def test_center_removes_the_mean_from_both_x_and_y():
    rng = np.random.default_rng(0)
    X = rng.normal(loc=5.0, scale=2.0, size=(50, 3))
    y = rng.normal(loc=-3.0, size=50)
    Xc, yc, X_mean, y_mean = r.center(X, y)
    assert np.allclose(Xc.mean(axis=0), 0.0, atol=1e-10)
    assert abs(float(yc.mean())) < 1e-10
    assert np.allclose(X_mean, X.mean(axis=0))
    assert y_mean == y.mean()


def test_make_dramatic_collinear_dataset_shape_and_true_coefficients():
    X, y, true_coef = r.make_dramatic_collinear_dataset(n=100, seed=0)
    assert X.shape == (100, 4)
    assert y.shape == (100,)
    assert np.array_equal(true_coef, np.array([1.0, 2.0, 3.0, 4.0]))
    # the fourth column is almost the first
    assert np.max(np.abs(X[:, 3] - X[:, 0])) < 1e-5


def test_normal_equations_and_lstsq_agree_on_a_well_conditioned_toy_problem():
    rng = np.random.default_rng(1)
    X = rng.normal(size=(200, 4))
    true_coef = np.array([1.0, -2.0, 0.5, 3.0])
    y = X @ true_coef + rng.normal(scale=0.01, size=200)
    beta_ne = r.fit_normal_equations(X, y)
    beta_lstsq = r.fit_lstsq(X, y)
    assert np.max(np.abs(beta_ne - beta_lstsq)) < 1e-8
    assert np.max(np.abs(beta_ne - true_coef)) < 0.01


def test_ols_regressor_exposes_sklearn_style_learned_attributes():
    rng = np.random.default_rng(2)
    X = rng.normal(size=(60, 3))
    y = X @ np.array([1.0, 2.0, -1.0]) + 4.0
    model = r.OLSRegressor(method="normal").fit(X, y)
    assert hasattr(model, "coef_")
    assert hasattr(model, "intercept_")
    assert hasattr(model, "n_features_in_")
    assert model.n_features_in_ == 3
    # get_params/set_params come from BaseEstimator's __init__ introspection
    params = model.get_params()
    assert params["method"] == "normal"
    assert params["fit_intercept"] is True

lesson_id: D153
day: 153
kind: guided-build
languages:
  - python
  - bash
setup_commands:
  - cd labs/sections/machine-learning/day-153-linear-regression-from-scratch
  - python3 -m venv .venv
  - .venv/bin/pip install -r requirements/requirements.txt
  - >-
    .venv/bin/python3 -c "import numpy, sklearn; print(numpy.__version__,
    sklearn.__version__)"
run_commands:
  - .venv/bin/pytest examples -q
  - .venv/bin/pytest starter -q
  - .venv/bin/python3 examples/report_measurements.py
test_commands:
  - bash tests/run_tests.sh
cleanup_commands:
  - >-
    find . -path ./.venv -prune -o -type d -name '__pycache__' -print -exec rm -rf -- {}
    +
  - rm -rf .pytest_cache
  - 'rm -rf .venv  # optional: removes the lab virtual environment'
  - 'git checkout -- starter/  # optional: reset your work'
requires_network: true
requires_api_key: false
estimated_minutes: 70
last_executed: '2026-08-27'
executed_on: >-
  macOS 26.5.2 (Apple Silicon, arm64, CPU only -- no GPU is needed or used), Python
  3.14.0, numpy 2.5.2, scikit-learn 1.9.0, pytest 9.1.1, bash 3.2.57 -- bash
  tests/run_tests.sh -> 14 checks, 0 failure(s), exit 0. pytest examples -q -> 15 passed.
  pytest starter -q -> 5 passed, 10 skipped (the five machinery checks in
  test_regression_lib.py are solved in both directories; the ten 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 is scikit-learn's own bundled load_diabetes, which ships inside the installed
  package and is never downloaded; 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 15 passed, rewrites the assertion
  `assert beta_sk[1] == pytest.approx(2.5, abs=0.05)` to expect 99999.0 instead, confirms
  a non-zero exit naming the failing test (test_02_a_near_duplicate_column_makes_the_
  normal_equations_explode), and removes the scratch directory. MEASURED PAIRS, all
  captured verbatim in expected-output/measured-values.txt. (1) On load_diabetes(scaled=
  True) with an intercept column prepended: max |normal equations - sklearn| = 1.2153e-10,
  max |lstsq - sklearn| = 1.1990e-12 -- lstsq about 101x closer to sklearn's own answer on
  the SAME well-conditioned data. (1b) condition number of X (with intercept) = 227.2248,
  condition number of X'X = 51631.1119, and cond(X'X)/cond(X)^2 = 1.0000000000, confirming
  the squaring relationship to ten decimal places. (2) On 100 rows, 3 random columns plus a
  fourth equal to column 0 plus 1e-7-scale noise, true coefficients [1,2,3,4], seed 0:
  normal-equation coefficients [0.001, 196747.976, 1.997, 2.994, -196742.975], lstsq
  coefficients [0.001, 207112.776, 1.997, 2.994, -207107.775], sklearn coefficients [0.001,
  2.501, 2.0, 2.997, 2.501] -- both closed-form routes explode past 196000 in opposite
  directions on the duplicated pair, sklearn's SVD-based minimum-norm solve splits the
  shared weight of 1 into 2.5 and 2.5 and recovers the other two coefficients near 2.0 and
  3.0. (2b) condition number of the same design matrix = 2.4363e+07, condition number of
  X'X = 5.6547e+14, and the measured ratio to cond(X)^2 is 0.9527, not 1.0 -- an honesty
  call, below. (3) On standardized diabetes features, the Day 111 stability threshold
  |1 - eta*a| < 1 evaluates to eta < 0.2485; at eta=0.2 (80 percent of threshold), gradient
  descent needs 3263 iterations to agree with the closed form to 3 decimals, 5277 for 6,
  and 7291 for 9 -- roughly linear growth, not exponential. The Hessian eigenvalue ratio
  (max/min) is 470.08 standardized. (3b) On RAW (scaled=False), centred diabetes features,
  the Hessian eigenvalue ratio is 76278.96 -- over 160x worse than standardized -- and the
  stability threshold shrinks to 4.8746e-04; at 95 percent of THAT threshold, after 200000
  iterations gradient descent has NOT reached even one decimal of agreement with the closed
  form (remaining max |coef - closed form| = 0.4692). (4) At 80 percent of the standardized
  threshold, gradient descent converges to 1e-9 in exactly 7132 iterations; at 102 percent
  of the same threshold, it diverges to non-finite values within 20000 iterations -- the
  Day 111 formula predicts exactly where that line falls. (5) Operation counts, no timing:
  the closed form on the diabetes shape (n=442, p=11 with intercept) costs n*p^2 + p^3 =
  54813 multiply-adds; gradient descent for the 7291 iterations exercise 3 measured for
  9-decimal agreement costs 2*n*p*iterations = 64452440 -- about 1176x more operations for
  the iterative method to match what one solve achieves. (6) OLSRegressor (inheriting
  BaseEstimator and RegressorMixin, per Day 146's measured requirement) matches sklearn's
  LinearRegression to within 4.4e-11 across all three fitting methods on standardized
  diabetes data. sklearn.utils.estimator_checks.check_estimator reports 52 total checks: 48
  passed, 2 failed by name (check_n_features_in_after_fitting, check_dtype_object -- both
  about input validation this implementation does not perform), 2 skipped by name
  (check_array_api_input: SCIPY_ARRAY_API not set; check_regressor_data_not_an_array:
  pandas not installed -- both environment conditions, not code defects). (7) fit_intercept
  by centring X and y and recovering the intercept afterwards, versus appending a column of
  ones and solving directly, agree to within 1.9554e-11 on the coefficients and 2.8422e-14
  on the intercept. TWO HONESTY CALLS. FIRST: on the dramatic near-duplicate-column case,
  the theoretical relationship cond(X'X) = cond(X)^2 -- verified to ten decimal places on
  the well-conditioned diabetes data in (1b) -- measures at only 0.9527 of that prediction
  here, because at cond(X) above ten million, computing the matrix's own smallest singular
  value is itself numerically imprecise; even the verification of the squaring relationship
  degrades under extreme enough ill-conditioning, and the lab reports the degraded ratio
  rather than asserting the clean theorem holds everywhere. SECOND: an early draft measured
  whether lstsq is closer to sklearn than the normal equations on ARBITRARY well-conditioned
  random data (cond(X) near 1.3) and found the ordering was a coin flip -- both routes sit
  at machine-epsilon noise when there is essentially nothing for a squared condition number
  to amplify. Harness check 8 therefore verifies the lstsq advantage on moderately
  ill-conditioned constructions (cond(X) in the 15-40 range, comparable to diabetes'
  227) across five seeds, where the effect is measured to hold reliably, rather than
  claiming it holds at any condition number whatsoever.

# Requirements

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

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

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

## Why the versions are pinned exactly

Every measured pair in this lab is deterministic given a seed -- there is
no repeated random sampling to average over, unlike Day 144's lab -- but
`np.linalg.lstsq`, `np.linalg.solve`, `np.linalg.eigvalsh` and
scikit-learn's own `LinearRegression` all bottom out in LAPACK routines
whose exact floating-point results can differ, in the last few bits, across
NumPy, SciPy and BLAS builds. Nothing here changes direction or order of
magnitude across ordinary machines, but the trailing decimals of numbers
like `1.2153e-10` are a property of this exact software stack, not of
mathematics.

What does not depend on the pins: every structural claim -- lstsq is closer
to sklearn than the normal equations are, `cond(X'X)` equals `cond(X)`
squared, the near-duplicate column makes the normal equations and lstsq
explode while sklearn stays sane, gradient descent diverges above the Day
111 stability threshold and converges below it, and the closed form uses
far fewer operations than gradient descent. `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 this lab uses is scikit-learn's bundled `load_diabetes`, which
ships inside the scikit-learn package itself and is never downloaded.

## Free and open-source status

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

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

# Day 153 lab brief -- Linear Regression from Scratch

Day 148 through 152 taught you what a fitted line means, where the formula
for its coefficients comes from, what happens with many predictors, and how
to regularize and score the result. Nobody in that run asked what
`LinearRegression().fit()` actually *does* when you call it.

This lab builds it three ways and finds out.

## The claim you are here to measure

> The textbook formula for linear regression is not what a good library
> actually runs, and you can show exactly why.

Exercise 1 fits ordinary least squares two ways -- the normal equations
`(X'X)^-1 X'y`, straight off the textbook page, and an `lstsq`-based solve
-- and checks both against scikit-learn's own `LinearRegression` on the
bundled diabetes dataset:

```text
max |normal equations - sklearn| : 1.2153e-10
max |lstsq            - sklearn| : 1.1990e-12
```

`lstsq` lands about a hundred times closer. Exercise 1b explains why:
`cond(X'X)` is exactly the square of `cond(X)` -- 51631.11 against 227.22
here -- and squaring a condition number is precisely how many digits of
precision the normal equations throw away that a direct solve of `X` does
not.

## The part that should genuinely surprise you

Exercise 2 makes the effect impossible to miss. Take three random columns,
add a fourth that is column 0 plus a sliver of noise -- 1e-7 in scale, far
below any real measurement error -- and fit with true coefficients
`[1, 2, 3, 4]`:

```text
normal-equation coefficients : [0.001, 196747.976, 1.997, 2.994, -196742.975]
lstsq coefficients           : [0.001, 207112.776, 1.997, 2.994, -207107.775]
sklearn coefficients         : [0.001, 2.501, 2.0, 2.997, 2.501]
```

Both from-scratch methods explode to plus and minus two hundred thousand.
sklearn stays sane, splitting the true value of 1 (which the near-duplicate
pair shares between them, since column 3 is almost column 0) into 2.5 and
2.5 -- an SVD-based minimum-norm solve that neither of the two textbook
routes performs.

## Gradient descent, and where Day 111 arrives

Exercise 3 fits the same problem with gradient descent instead, on
standardized features, and measures exactly how many iterations it takes
to agree with the closed form to 3, 6 and 9 decimal places. Exercise 3b
repeats it on the RAW, unscaled diabetes columns -- age in years next to
sex coded 1 or 2 next to serum measurements in the hundreds -- and shows
gradient descent barely moving, because Day 111's condition number (the
ratio of the loss's largest to smallest Hessian eigenvalue) is over a
hundred times worse unscaled. Exercise 4 finds the exact learning rate
where gradient descent stops converging and starts diverging, and shows it
lines up with Day 111's stability formula `|1 - eta * a| < 1` to the
decimal.

## What each method costs, without a stopwatch

Exercise 5 counts operations instead of timing anything -- the closed form
needs `n*p^2 + p^3` multiply-adds; gradient descent needs `2*n*p` per
iteration. On the diabetes shape, the closed form wins by a factor of over
a thousand, for the number of iterations exercise 3 measured.

## The estimator, and the two things it does not do

Exercise 6 wraps all three fitting routes in an `OLSRegressor` that
inherits `BaseEstimator` and `RegressorMixin` -- Day 146 already measured
why a from-scratch estimator needs that inheritance to survive
`Pipeline.predict()` and `cross_val_score`, and this lesson does not
re-teach it. Then it runs `sklearn.utils.estimator_checks.check_estimator`
against it and reports the real result: 48 of 52 checks pass, and two fail
by name, both about input validation this implementation does not perform.

## 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 ten 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.

"""Ordinary least squares, three ways: what the library was doing for you.

Three implementations of the same fit -- the normal equations, an
lstsq-based solve, and gradient descent -- plus a scikit-learn-compatible
estimator wrapping all three, and then a direct measurement of where they
agree and where they do not.

Day 111 already derived gradient descent, its update rule, the stability
condition ``|1 - eta * a| < 1`` for an eigenvalue ``a`` of the loss's
Hessian, and the condition number as the ratio of the Hessian's largest and
smallest eigenvalues. None of that is re-derived here; it is applied.

Day 146 already established the scikit-learn estimator API contract --
``fit``/``predict``, learned attributes with a trailing underscore -- and
measured that a from-scratch estimator raises ``AttributeError`` on
``__sklearn_tags__`` inside ``Pipeline.predict()`` and ``cross_val_score``
unless it inherits ``BaseEstimator``. ``OLSRegressor`` below inherits it for
exactly that reason, without re-teaching the contract.

Everything here is deterministic given a seed.
"""

from __future__ import annotations

import numpy as np

from sklearn.base import BaseEstimator, RegressorMixin
from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression
from sklearn.utils.validation import check_array, check_is_fitted, check_X_y

# --------------------------------------------------------------------------
# 1. Loading the data, and the two condition numbers
# --------------------------------------------------------------------------


def load_diabetes_data(scaled: bool = True):
    """The bundled diabetes regression dataset: X (442, 10), y (442,).

    ``scaled=True`` (the default) gives sklearn's own mean-centred,
    unit-L2-norm-scaled columns. ``scaled=False`` gives the raw clinical
    units -- age in years, sex coded 1/2, six serum measurements on very
    different scales -- which is the badly scaled version used to show
    gradient descent struggling.
    """
    data = load_diabetes(scaled=scaled)
    return data.data, data.target


def add_intercept_column(X: np.ndarray) -> np.ndarray:
    """Prepend a column of ones, the classic way to fold in an intercept."""
    ones = np.ones((X.shape[0], 1))
    return np.hstack([ones, X])


def condition_numbers(A: np.ndarray) -> tuple[float, float]:
    """``(cond(A), cond(A'A))`` -- the second is exactly the square of the first.

    ``cond`` here is the ratio of largest to smallest singular value. Squaring
    a matrix's condition number is the textbook reason the normal equations
    lose precision that a direct solve of ``A`` does not.
    """
    cond_a = float(np.linalg.cond(A))
    cond_ata = float(np.linalg.cond(A.T @ A))
    return cond_a, cond_ata


# --------------------------------------------------------------------------
# 2. Three ways to fit: normal equations, lstsq, and sklearn as the referee
# --------------------------------------------------------------------------


def fit_normal_equations(A: np.ndarray, y: np.ndarray) -> np.ndarray:
    """Solve the normal equations ``A'A b = A'y`` directly.

    This is the formula from the textbook page: form ``A'A``, form ``A'y``,
    solve the square system. It is also, as this module measures, the least
    numerically careful of the three routes here.
    """
    return np.linalg.solve(A.T @ A, A.T @ y)


def fit_lstsq(A: np.ndarray, y: np.ndarray) -> np.ndarray:
    """Solve the same least-squares problem via ``numpy.linalg.lstsq``.

    ``lstsq`` factors ``A`` directly (an SVD-based solve under the hood)
    rather than forming and inverting ``A'A``, which is why it does not
    inherit the squared condition number of the normal equations.
    """
    coef, *_ = np.linalg.lstsq(A, y, rcond=None)
    return coef


def sklearn_reference_fit(X: np.ndarray, y: np.ndarray) -> np.ndarray:
    """Fit scikit-learn's own ``LinearRegression`` and return ``[intercept, *coef]``.

    Used as the referee throughout this module: not because it is a fourth
    algorithm, but because it is the thing every from-scratch attempt here is
    being measured against.
    """
    model = LinearRegression().fit(X, y)
    return np.concatenate([[float(model.intercept_)], model.coef_])


def max_abs_difference(a: np.ndarray, b: np.ndarray) -> float:
    return float(np.max(np.abs(a - b)))


# --------------------------------------------------------------------------
# 3. The dramatic case: a near-duplicate column
# --------------------------------------------------------------------------


def make_dramatic_collinear_dataset(n: int = 100, seed: int = 0, noise_scale: float = 1e-7):
    """Three random predictors plus a fourth that is almost the first.

    The fourth column is column 0 plus a sliver of noise -- 1e-7 in scale,
    far smaller than any real measurement error. The true coefficients are
    [1, 2, 3, 4]; a well-conditioned fit should recover something close to
    that. It does not, for two of the three methods below.
    """
    rng = np.random.default_rng(seed)
    X = rng.normal(size=(n, 3))
    duplicate = X[:, 0] + rng.normal(scale=noise_scale, size=n)
    X = np.column_stack([X, duplicate])
    true_coef = np.array([1.0, 2.0, 3.0, 4.0])
    y = X @ true_coef + rng.normal(scale=0.1, size=n)
    return X, y, true_coef


# --------------------------------------------------------------------------
# 4. Gradient descent, and the Day 111 stability threshold arriving here
# --------------------------------------------------------------------------


def standardize(X: np.ndarray) -> np.ndarray:
    """Zero mean, unit standard deviation, column by column."""
    mu = X.mean(axis=0)
    sd = X.std(axis=0)
    return (X - mu) / sd


def center(X: np.ndarray, y: np.ndarray):
    """Subtract each column's mean from X and y's mean from y.

    Centring removes the need for a separate intercept term: the fitted
    hyperplane through centred data passes through the origin, and the
    intercept is recovered afterwards as ``y.mean() - X.mean(axis=0) @ coef``.
    """
    X_mean = X.mean(axis=0)
    y_mean = y.mean()
    return X - X_mean, y - y_mean, X_mean, y_mean


def hessian_eigenvalues(X: np.ndarray, n: int) -> np.ndarray:
    """Eigenvalues of the mean-squared-error loss's Hessian, ``(2/n) X'X``.

    Day 111 established that gradient descent on a quadratic loss is stable
    exactly when ``|1 - eta * a| < 1`` for every eigenvalue ``a`` of this
    Hessian, and that the ratio of its largest to smallest eigenvalue is the
    condition number governing how slowly the slowest direction converges.
    """
    XtX = X.T @ X
    return np.linalg.eigvalsh((2.0 / n) * XtX)


def stability_threshold(X: np.ndarray) -> float:
    """The largest learning rate for which every eigenvalue keeps |1 - eta*a| < 1.

    From Day 111's condition: stability requires ``eta < 2 / a`` for every
    eigenvalue ``a``, so the binding constraint is the largest eigenvalue.
    """
    n = X.shape[0]
    eig_max = float(hessian_eigenvalues(X, n).max())
    return 2.0 / eig_max


def fit_gradient_descent(X: np.ndarray, y: np.ndarray, lr: float, n_iter: int) -> np.ndarray:
    """Plain batch gradient descent on the centred mean-squared-error loss.

    Starts at the zero vector. ``X`` and ``y`` are assumed already centred,
    so no intercept term appears in this loop -- see ``center`` above.
    """
    n, p = X.shape
    coef = np.zeros(p)
    for _ in range(n_iter):
        grad = (2.0 / n) * (X.T @ (X @ coef - y))
        coef = coef - lr * grad
    return coef


def iters_to_tolerance(X: np.ndarray, y: np.ndarray, lr: float, target: np.ndarray, tol: float, max_iter: int):
    """How many gradient-descent iterations until every coefficient is within tol of target.

    Returns ``(iterations, final_coef)``, or ``(None, final_coef)`` if the
    tolerance was never reached within ``max_iter`` -- either because
    convergence needs more steps, or because the run diverged. A run that
    stops producing finite numbers returns ``("diverged", coef)``.
    """
    n, p = X.shape
    coef = np.zeros(p)
    for i in range(1, max_iter + 1):
        grad = (2.0 / n) * (X.T @ (X @ coef - y))
        coef = coef - lr * grad
        if not np.all(np.isfinite(coef)):
            return "diverged", coef
        if float(np.max(np.abs(coef - target))) < tol:
            return i, coef
    return None, coef


# --------------------------------------------------------------------------
# 5. Counting operations instead of timing them
# --------------------------------------------------------------------------


def normal_equation_op_count(n: int, p: int) -> int:
    """Multiply-adds to form A'A (n*p^2) plus solve the p-by-p system (p^3)."""
    return n * p * p + p**3


def gradient_descent_op_count(n: int, p: int, iterations: int) -> int:
    """Multiply-adds for `iterations` steps: two n*p matrix-vector products each."""
    return 2 * n * p * iterations


# --------------------------------------------------------------------------
# 6. A scikit-learn-compatible estimator
# --------------------------------------------------------------------------


class OLSRegressor(RegressorMixin, BaseEstimator):
    """Ordinary least squares, fit by one of three methods, as a real estimator.

    Inherits ``RegressorMixin`` and ``BaseEstimator`` because Day 146 already
    measured what happens without them: ``fit``, ``predict`` and ``score``
    work perfectly well called directly, but ``Pipeline.predict()`` and
    ``cross_val_score`` both raise ``AttributeError`` on ``__sklearn_tags__``,
    which only ``BaseEstimator`` supplies. That lesson is assumed, not
    repeated.

    ``fit_intercept=True`` centres X and y, fits the centred problem, and
    recovers the intercept afterwards -- rather than appending a column of
    ones -- because the two are measured elsewhere in this module to agree
    to about ten decimal places and centring avoids growing the design
    matrix by one column.
    """

    def __init__(self, method: str = "lstsq", fit_intercept: bool = True, lr: float = 0.1, n_iter: int = 1000):
        self.method = method
        self.fit_intercept = fit_intercept
        self.lr = lr
        self.n_iter = n_iter

    def fit(self, X, y):
        X, y = check_X_y(X, y)
        if self.fit_intercept:
            Xc, yc, X_mean, y_mean = center(X, y)
        else:
            Xc, yc = X, y
            X_mean, y_mean = np.zeros(X.shape[1]), 0.0

        if self.method == "normal":
            coef = fit_normal_equations(Xc, yc)
        elif self.method == "lstsq":
            coef = fit_lstsq(Xc, yc)
        elif self.method == "gd":
            coef = fit_gradient_descent(Xc, yc, self.lr, self.n_iter)
        else:
            raise ValueError(f"unknown method {self.method!r}; use 'normal', 'lstsq' or 'gd'")

        self.coef_ = coef
        self.intercept_ = float(y_mean - X_mean @ coef) if self.fit_intercept else 0.0
        self.n_features_in_ = X.shape[1]
        return self

    def predict(self, X):
        check_is_fitted(self, "coef_")
        X = check_array(X)
        return X @ self.coef_ + self.intercept_


def run_check_estimator(estimator):
    """Run scikit-learn's estimator_checks suite and tally results by name.

    Returns ``(passed, failed, skipped)`` where ``failed`` and ``skipped``
    are lists of ``(check_name, message)`` pairs -- never silently discarded.
    """
    from sklearn.utils.estimator_checks import check_estimator

    results = []

    def record(*, estimator, check_name, exception, status, expected_to_fail, expected_to_fail_reason):
        results.append((check_name, status, str(exception)[:200] if exception else None))

    check_estimator(estimator, on_fail=None, on_skip=None, callback=record)
    passed = [name for name, status, _msg in results if status == "passed"]
    failed = [(name, msg) for name, status, msg in results if status == "failed"]
    skipped = [(name, msg) for name, status, msg in results if status == "skipped"]
    return passed, failed, skipped


# --------------------------------------------------------------------------
# 7. fit_intercept two ways: centring versus an appended column
# --------------------------------------------------------------------------


def fit_intercept_two_ways(X: np.ndarray, y: np.ndarray):
    """Compare centring against appending a ones column, on the same data.

    Returns ``(coef_column, intercept_column, coef_centred, intercept_centred)``.
    """
    n = X.shape[0]
    A = add_intercept_column(X)
    beta_column = fit_normal_equations(A, y)
    intercept_column, coef_column = float(beta_column[0]), beta_column[1:]

    Xc, yc, X_mean, y_mean = center(X, y)
    coef_centred = fit_normal_equations(Xc, yc)
    intercept_centred = float(y_mean - X_mean @ coef_centred)

    return coef_column, intercept_column, coef_centred, intercept_centred

"""Ten exercises in what a linear-regression library was doing for you.

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_diabetes_data(scaled=True)


@pytest.fixture(scope="module")
def diabetes_raw():
    return r.load_diabetes_data(scaled=False)


def test_01_lstsq_agrees_with_sklearn_a_hundred_times_more_closely_than_the_normal_equations(diabetes):
    pytest.skip(
        "Prepend an intercept column with r.add_intercept_column, fit "
        "r.fit_normal_equations and r.fit_lstsq on it, and fit "
        "r.sklearn_reference_fit on the raw X. Assert "
        "r.max_abs_difference(normal_equations, sklearn) is about 1.2153e-10 "
        "and the lstsq gap is about 1.199e-12 -- roughly a hundred times "
        "closer to sklearn's own answer, on the SAME well-conditioned data."
    )


def test_01b_the_normal_equations_condition_number_is_exactly_the_square(diabetes):
    pytest.skip(
        "Call r.condition_numbers on the intercept-augmented design matrix. "
        "Assert cond(A) is about 227.2248 and cond(A'A) is about 51631.1119, "
        "and that cond(A'A) / cond(A)**2 is within 1e-8 of 1.0. This is the "
        "textbook reason the normal equations lose precision that a direct "
        "solve of A does not."
    )


def test_02_a_near_duplicate_column_makes_the_normal_equations_explode():
    pytest.skip(
        "Build r.make_dramatic_collinear_dataset(n=100, seed=0): three "
        "random columns plus a fourth that is column 0 plus a sliver of "
        "noise, true coefficients [1, 2, 3, 4]. Fit both closed forms and "
        "sklearn's own LinearRegression. Assert the normal-equation and "
        "lstsq coefficients for the duplicated pair explode past 1e5 in "
        "opposite directions, while sklearn's stay near 2.5 and 2.5 -- an "
        "SVD-based minimum-norm solve splitting the weight evenly -- and "
        "sklearn recovers the other two coefficients near 2.0 and 3.0."
    )


def test_02b_even_the_squaring_relationship_becomes_hard_to_verify_here():
    pytest.skip(
        "On the same dramatic dataset, assert cond(A) exceeds 1e6 and "
        "cond(A'A) exceeds 1e13 -- several orders of magnitude worse than "
        "the diabetes case. Then assert cond(A'A) / cond(A)**2 lies strictly "
        "between 0.9 and 1.0 but differs from 1.0 by MORE than 1e-4: at this "
        "level of ill-conditioning, computing the smallest singular value of "
        "an already near-singular matrix is itself imprecise, so even the "
        "squaring relationship's own verification degrades."
    )


def test_03_gradient_descent_reaches_the_closed_form_at_three_six_and_nine_decimals(diabetes):
    pytest.skip(
        "Standardize X with r.standardize, centre y, and fit the closed "
        "form as the target. Compute r.stability_threshold and assert it is "
        "about 0.2485 (Day 111's condition applied to this Hessian). At a "
        "learning rate of 0.2 -- about 80 percent of that threshold -- "
        "assert r.iters_to_tolerance needs exactly 3263 iterations for 3 "
        "decimals, 5277 for 6, and 7291 for 9. Then assert the growth is "
        "roughly linear, not exponential: iters_6 - iters_3 < iters_3."
    )


def test_03b_the_same_setup_on_unscaled_features_barely_moves(diabetes, diabetes_raw):
    pytest.skip(
        "Compute r.hessian_eigenvalues for the standardized diabetes data "
        "and for the RAW (scaled=False), centred diabetes data. Assert the "
        "max/min eigenvalue ratio is about 470.08 standardized and about "
        "76278.96 raw -- over a hundred times worse conditioned, in exactly "
        "the sense Day 111 defined. Then, at 95 percent of the raw data's "
        "own stability threshold, assert r.iters_to_tolerance has NOT "
        "converged after 200,000 iterations and the remaining gap from the "
        "closed form still exceeds 0.1 -- a learning rate that is stable in "
        "principle can still be catastrophically slow."
    )


def test_04_the_day_111_stability_threshold_predicts_divergence_exactly(diabetes):
    pytest.skip(
        "At 80 percent of r.stability_threshold, assert "
        "r.iters_to_tolerance converges (to 1e-9) in exactly 7132 "
        "iterations. At 102 percent of the same threshold, assert it "
        "returns the string 'diverged' and the returned coefficients are "
        "no longer all finite. The formula from Day 111, "
        "|1 - eta * a| < 1, predicts exactly where that line falls."
    )


def test_05_the_closed_form_needs_a_thousand_times_fewer_operations(diabetes):
    pytest.skip(
        "Using the diabetes shape (n=442, p=10), compute "
        "r.normal_equation_op_count(442, 11) -- the +1 is the intercept "
        "column -- and assert it equals 54813. Compute "
        "r.gradient_descent_op_count(442, 10, 7291) using exercise 3's "
        "9-decimal iteration count and assert it equals 64452440. Assert "
        "the ratio is between 1000 and 1300: over a thousand times more "
        "multiply-adds for gradient descent to match what the closed form "
        "gets in one solve. No wall-clock timing anywhere -- count "
        "operations instead."
    )


def test_06_the_estimator_matches_sklearn_and_check_estimator_names_two_failures(diabetes):
    pytest.skip(
        "Fit r.OLSRegressor with method='normal', 'lstsq' and 'gd' "
        "(lr=0.2, n_iter=8000) on standardized diabetes data, and assert "
        "each one's coef_ matches sklearn's own LinearRegression to within "
        "1e-8 or better. Then call r.run_check_estimator(r.OLSRegressor()) "
        "and assert exactly 48 of 52 checks pass, that the two failures are "
        "named check_n_features_in_after_fitting and check_dtype_object, "
        "and that the two skips are check_array_api_input and "
        "check_regressor_data_not_an_array -- do not suppress or hide "
        "either list."
    )


def test_07_centring_and_appending_a_column_agree_to_ten_decimal_places(diabetes):
    pytest.skip(
        "Call r.fit_intercept_two_ways on the diabetes data and assert the "
        "max coefficient difference between the column-append approach and "
        "the centring approach is below 1e-9, and the intercept difference "
        "is below 1e-9 as well. Two routes to the same intercept, agreeing "
        "to nine decimal places or better."
    )

"""Machinery checks -- already solved, in both starter/ and examples/.

These confirm the library behaves as documented. They are not the
exercises; `test_regression_claims.py` is.
"""

import numpy as np

import regression_lib as r


def test_add_intercept_column_prepends_a_column_of_ones():
    X = np.arange(6.0).reshape(3, 2)
    A = r.add_intercept_column(X)
    assert A.shape == (3, 3)
    assert np.array_equal(A[:, 0], np.ones(3))
    assert np.array_equal(A[:, 1:], X)


def test_center_removes_the_mean_from_both_x_and_y():
    rng = np.random.default_rng(0)
    X = rng.normal(loc=5.0, scale=2.0, size=(50, 3))
    y = rng.normal(loc=-3.0, size=50)
    Xc, yc, X_mean, y_mean = r.center(X, y)
    assert np.allclose(Xc.mean(axis=0), 0.0, atol=1e-10)
    assert abs(float(yc.mean())) < 1e-10
    assert np.allclose(X_mean, X.mean(axis=0))
    assert y_mean == y.mean()


def test_make_dramatic_collinear_dataset_shape_and_true_coefficients():
    X, y, true_coef = r.make_dramatic_collinear_dataset(n=100, seed=0)
    assert X.shape == (100, 4)
    assert y.shape == (100,)
    assert np.array_equal(true_coef, np.array([1.0, 2.0, 3.0, 4.0]))
    # the fourth column is almost the first
    assert np.max(np.abs(X[:, 3] - X[:, 0])) < 1e-5


def test_normal_equations_and_lstsq_agree_on_a_well_conditioned_toy_problem():
    rng = np.random.default_rng(1)
    X = rng.normal(size=(200, 4))
    true_coef = np.array([1.0, -2.0, 0.5, 3.0])
    y = X @ true_coef + rng.normal(scale=0.01, size=200)
    beta_ne = r.fit_normal_equations(X, y)
    beta_lstsq = r.fit_lstsq(X, y)
    assert np.max(np.abs(beta_ne - beta_lstsq)) < 1e-8
    assert np.max(np.abs(beta_ne - true_coef)) < 0.01


def test_ols_regressor_exposes_sklearn_style_learned_attributes():
    rng = np.random.default_rng(2)
    X = rng.normal(size=(60, 3))
    y = X @ np.array([1.0, 2.0, -1.0]) + 4.0
    model = r.OLSRegressor(method="normal").fit(X, y)
    assert hasattr(model, "coef_")
    assert hasattr(model, "intercept_")
    assert hasattr(model, "n_features_in_")
    assert model.n_features_in_ == 3
    # get_params/set_params come from BaseEstimator's __init__ introspection
    params = model.get_params()
    assert params["method"] == "normal"
    assert params["fit_intercept"] is True

#!/usr/bin/env bash
# Day 153 lab harness: "Linear Regression from Scratch"
#
# 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
from sklearn.linear_model import LinearRegression

import regression_lib as r

errors = []


def expect(label, got, want, tol=1e-6):
    if isinstance(want, float):
        if abs(got - want) > tol:
            errors.append(f"{label}: expected {want}, got {got}")
    elif got != want:
        errors.append(f"{label}: expected {want}, got {got}")


X, y = r.load_diabetes_data(scaled=True)
n, p = X.shape

# 1. Three closed forms on well-conditioned data
A = r.add_intercept_column(X)
beta_ne = r.fit_normal_equations(A, y)
beta_lstsq = r.fit_lstsq(A, y)
beta_sk = r.sklearn_reference_fit(X, y)
gap_ne = r.max_abs_difference(beta_ne, beta_sk)
gap_lstsq = r.max_abs_difference(beta_lstsq, beta_sk)
if not (gap_ne / gap_lstsq > 50):
    errors.append(f"lstsq was not roughly a hundred times closer: ratio {gap_ne / gap_lstsq}")

# 1b. cond(X'X) is the square of cond(X)
cond_a, cond_ata = r.condition_numbers(A)
expect("cond(X) with intercept", round(cond_a, 4), 227.2248, tol=0.01)
expect("cond(X'X)", round(cond_ata, 4), 51631.1119, tol=0.5)
if abs(cond_ata / cond_a**2 - 1.0) > 1e-6:
    errors.append(f"cond(X'X)/cond(X)^2 was not ~1.0: {cond_ata / cond_a**2}")

# 2. The near-duplicate column
Xd, yd, true_coef = r.make_dramatic_collinear_dataset(n=100, seed=0)
Ad = r.add_intercept_column(Xd)
beta_ne_d = r.fit_normal_equations(Ad, yd)
beta_lstsq_d = r.fit_lstsq(Ad, yd)
beta_sk_d = r.sklearn_reference_fit(Xd, yd)
if not (beta_ne_d[1] > 1e5 and beta_ne_d[4] < -1e5):
    errors.append("normal equations did not explode on the near-duplicate column")
if not (beta_lstsq_d[1] > 1e5 and beta_lstsq_d[4] < -1e5):
    errors.append("lstsq did not explode on the near-duplicate column")
if not (abs(beta_sk_d[1] - 2.5) < 0.05 and abs(beta_sk_d[4] - 2.5) < 0.05):
    errors.append("sklearn did not split the weight evenly near 2.5 and 2.5")

# 2b. Squaring relationship degrades under extreme ill-conditioning
cond_ad, cond_atad = r.condition_numbers(Ad)
if not (cond_ad > 1e6 and cond_atad > 1e13):
    errors.append("dramatic-case condition numbers were not as extreme as expected")
ratio_d = cond_atad / cond_ad**2
if not (0.9 < ratio_d < 1.0 and abs(ratio_d - 1.0) > 1e-4):
    errors.append(f"squaring relationship did not degrade as expected: ratio {ratio_d}")

# 3. Gradient descent vs closed form, standardized
Xs = r.standardize(X)
yc = y - y.mean()
target = r.fit_normal_equations(Xs, yc)
threshold = r.stability_threshold(Xs)
expect("stability threshold, standardized", round(threshold, 4), 0.2485, tol=0.001)
iters_3, _ = r.iters_to_tolerance(Xs, yc, 0.2, target, 5e-4, 200_000)
iters_6, _ = r.iters_to_tolerance(Xs, yc, 0.2, target, 5e-7, 200_000)
iters_9, _ = r.iters_to_tolerance(Xs, yc, 0.2, target, 5e-10, 200_000)
expect("iterations for 3 decimals", iters_3, 3263)
expect("iterations for 6 decimals", iters_6, 5277)
expect("iterations for 9 decimals", iters_9, 7291)

# 3b. Raw features
Xraw, yraw = r.load_diabetes_data(scaled=False)
Xrc, yrc, _, _ = r.center(Xraw, yraw)
eig_s = r.hessian_eigenvalues(Xs, n)
eig_r = r.hessian_eigenvalues(Xrc, n)
ratio_scaled = float(eig_s.max() / eig_s.min())
ratio_raw = float(eig_r.max() / eig_r.min())
expect("Hessian eigenvalue ratio, standardized", round(ratio_scaled, 2), 470.08, tol=0.5)
expect("Hessian eigenvalue ratio, raw", round(ratio_raw, 2), 76278.96, tol=5.0)
raw_threshold = r.stability_threshold(Xrc)
target_r = r.fit_normal_equations(Xrc, yrc)
status_r, coef_r = r.iters_to_tolerance(Xrc, yrc, raw_threshold * 0.95, target_r, 5e-4, 200_000)
if status_r is not None:
    errors.append("raw-feature gradient descent converged when it was expected to remain slow")
if r.max_abs_difference(coef_r, target_r) <= 0.1:
    errors.append("raw-feature gradient descent got closer than expected after 200000 iterations")

# 4. Stability threshold predicts divergence exactly
below_status, _ = r.iters_to_tolerance(Xs, yc, threshold * 0.8, target, 1e-9, 20_000)
above_status, above_coef = r.iters_to_tolerance(Xs, yc, threshold * 1.02, target, 1e-9, 20_000)
expect("iterations to converge at 80 percent of threshold", below_status, 7132)
if above_status != "diverged" or np.all(np.isfinite(above_coef)):
    errors.append("gradient descent did not diverge above the stability threshold")

# 5. Operation counts
ops_normal = r.normal_equation_op_count(n, p + 1)
ops_gd = r.gradient_descent_op_count(n, p, iters_9)
expect("normal-equation operation count", ops_normal, 54813)
expect("gradient-descent operation count", ops_gd, 64_452_440)

# 6. The estimator and check_estimator
sk = LinearRegression().fit(Xs, y)
normal_est = r.OLSRegressor(method="normal").fit(Xs, y)
lstsq_est = r.OLSRegressor(method="lstsq").fit(Xs, y)
gd_est = r.OLSRegressor(method="gd", lr=0.2, n_iter=8000).fit(Xs, y)
if r.max_abs_difference(normal_est.coef_, sk.coef_) > 1e-8:
    errors.append("OLSRegressor(method='normal') did not match sklearn closely enough")
if r.max_abs_difference(lstsq_est.coef_, sk.coef_) > 1e-8:
    errors.append("OLSRegressor(method='lstsq') did not match sklearn closely enough")
if r.max_abs_difference(gd_est.coef_, sk.coef_) > 1e-7:
    errors.append("OLSRegressor(method='gd') did not match sklearn closely enough")

passed, failed, skipped = r.run_check_estimator(r.OLSRegressor())
expect("check_estimator passed count", len(passed), 48)
failed_names = sorted(name for name, _msg in failed)
skipped_names = sorted(name for name, _msg in skipped)
expect("check_estimator failed names", failed_names, ["check_dtype_object", "check_n_features_in_after_fitting"])
expect(
    "check_estimator skipped names",
    skipped_names,
    ["check_array_api_input", "check_regressor_data_not_an_array"],
)

# 7. fit_intercept two ways
coef_col, intercept_col, coef_centred, intercept_centred = r.fit_intercept_two_ways(X, y)
if r.max_abs_difference(coef_col, coef_centred) > 1e-9:
    errors.append("centring and appending a column did not agree on coefficients")
if abs(intercept_col - intercept_centred) > 1e-9:
    errors.append("centring and appending a column did not agree on the intercept")

if errors:
    for e in errors:
        print("ERROR:", e)
    sys.exit(1)
print("all direct checks passed")
PYEOF
)
if echo "$DIRECT_CHECK" | grep -q "all direct checks passed"; then
  ok "exercises 1-7 reproduced directly against 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 "^15 passed"; then
  ok "pytest examples -q -> 15 passed"
else
  fail "pytest examples -q did not report 15 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, 10 skipped"; then
  ok "pytest starter -q -> 5 passed, 10 skipped (the machinery checks pass; the ten exercises are stubs)"
else
  fail "pytest starter -q did not report 5 passed, 10 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}/d153-scratch.XXXXXX")
cp examples/*.py "$SCRATCH"/
SCRATCH_OUT=$("$PYTEST" "$SCRATCH" -q 2>&1)
if echo "$SCRATCH_OUT" | tail -1 | grep -qE "^15 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 beta_sk[1] == pytest.approx(2.5, abs=0.05)"
replacement = "assert beta_sk[1] == pytest.approx(99999.0, abs=0.05)"
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_02_a_near_duplicate_column_makes_the_normal_equations_explode"; then
  ok "breaking exercise 2'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 shapes 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 = []

# The near-duplicate-column story is not a property of one seed.
for seed in range(5):
    X, y, true_coef = r.make_dramatic_collinear_dataset(n=100, seed=seed)
    A = r.add_intercept_column(X)
    beta_ne = r.fit_normal_equations(A, y)
    beta_sk = r.sklearn_reference_fit(X, y)
    if not (abs(beta_ne[1]) > 1e4 or abs(beta_ne[4]) > 1e4):
        problems.append(f"seed {seed}: normal equations did not explode on the duplicate pair")
    if r.max_abs_difference(beta_sk[[2, 3]], true_coef[[1, 2]]) > 0.2:
        problems.append(f"seed {seed}: sklearn did not recover the two clean coefficients")

# lstsq stays closer to sklearn than the normal equations, on OTHER
# moderately ill-conditioned datasets too -- not just the diabetes case.
# On near-orthogonal columns (cond close to 1) both routes sit at machine
# epsilon and the ordering is noise, which is itself a finding: the
# advantage specifically shows up once squaring the condition number
# starts to bite, not universally.
for seed in range(5):
    rng = np.random.default_rng(100 + seed)
    base = rng.normal(size=(300, 6))
    mix = np.full((6, 6), 0.995)
    np.fill_diagonal(mix, 1.0)
    X = base @ np.linalg.cholesky(mix).T
    true_coef = rng.normal(size=6)
    y = X @ true_coef + rng.normal(scale=0.05, size=300)
    A = r.add_intercept_column(X)
    beta_ne = r.fit_normal_equations(A, y)
    beta_lstsq = r.fit_lstsq(A, y)
    beta_sk = r.sklearn_reference_fit(X, y)
    gap_ne = r.max_abs_difference(beta_ne, beta_sk)
    gap_lstsq = r.max_abs_difference(beta_lstsq, beta_sk)
    if gap_lstsq > gap_ne:
        problems.append(f"seed {seed}: lstsq was not at least as close to sklearn as the normal equations")

# cond(X'X) equals cond(X) squared on well-conditioned matrices in general.
for seed in range(5):
    rng = np.random.default_rng(200 + seed)
    A = rng.normal(size=(150, 5))
    cond_a, cond_ata = r.condition_numbers(A)
    if abs(cond_ata / cond_a**2 - 1.0) > 1e-6:
        problems.append(f"seed {seed}: cond(A'A) was not the square of cond(A)")

if problems:
    for p in problems:
        print("ERROR:", p)
else:
    print("every direction held")
PYEOF
)
if [ "$DIRECTION" = "every direction held" ]; then
  ok "the duplicate-column explosion, lstsq's advantage, and the squaring relationship hold beyond the quoted seeds"
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