Hands-on lab — Day 149: Loss Functions and Least Squares

"""Loss functions, measured: what choosing squared error over absolute error
actually decides, and what each is implicitly betting on.

A loss function is a choice, not a law of nature. This module measures four
consequences of that choice: the minimiser of squared error is the mean and
the minimiser of absolute error is the median; the squared-error landscape is
smooth with one minimum while the absolute-error landscape is piecewise
linear and kinked; the normal equations solve squared error in closed form
because that smoothness gives a zero derivative one can solve for directly,
while absolute error has no such closed form; and swapping which loss you
minimise changes how far a single outlier can move your line, and which loss
wins depends on what the errors actually look like -- Gaussian, or
heavy-tailed.

Everything here is deterministic given a seed.
"""

from __future__ import annotations

import numpy as np

from sklearn.linear_model import HuberRegressor, LinearRegression, QuantileRegressor


# --------------------------------------------------------------------------
# 1. What each loss minimises
# --------------------------------------------------------------------------


def sse(values, candidate: float) -> float:
    """Sum of squared error between a scalar candidate and every value."""
    values = np.asarray(values, dtype=float)
    return float(np.sum((values - candidate) ** 2))


def sae(values, candidate: float) -> float:
    """Sum of absolute error between a scalar candidate and every value."""
    values = np.asarray(values, dtype=float)
    return float(np.sum(np.abs(values - candidate)))


def grid_minimize(values, loss_fn, lo: float, hi: float, steps: int = 200_001) -> float:
    """Find the candidate in a fine grid that minimises the given loss.

    A brute-force numerical stand-in for calculus: since it searches a
    finite grid it lands within ``(hi - lo) / (steps - 1)`` of the true
    minimiser, not exactly on it.
    """
    grid = np.linspace(lo, hi, steps)
    losses = np.array([loss_fn(values, c) for c in grid])
    return float(grid[int(np.argmin(losses))])


# --------------------------------------------------------------------------
# 2. The shape of the loss landscape: smooth, or kinked
# --------------------------------------------------------------------------


def make_line_data(n: int = 150, seed: int = 0, true_intercept: float = 5.0,
                    true_slope: float = 3.0, noise_sd: float = 2.0,
                    heavy_tailed: bool = False, heavy_df: int = 3,
                    heavy_scale: float = 1.2):
    """A simple one-predictor dataset with a known true line.

    ``heavy_tailed=True`` swaps the Gaussian error for a scaled Student's t
    with ``heavy_df`` degrees of freedom, which has the same rough central
    spread but far fatter tails -- the construction used to measure what
    squared error implicitly assumes about the errors.
    """
    rng = np.random.default_rng(seed)
    x = rng.uniform(0.0, 10.0, n)
    if heavy_tailed:
        errors = rng.standard_t(df=heavy_df, size=n) * heavy_scale
    else:
        errors = rng.normal(0.0, noise_sd, n)
    y = true_intercept + true_slope * x + errors
    return x, y


def loss_landscape(x, y, intercept: float, slopes):
    """Total squared error and total absolute error at each candidate slope,
    with the intercept held fixed. Returns ``(sq_losses, abs_losses)``.
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    slopes = np.asarray(slopes, dtype=float)
    sq_losses = np.empty_like(slopes)
    abs_losses = np.empty_like(slopes)
    for i, m in enumerate(slopes):
        residual = y - (intercept + m * x)
        sq_losses[i] = np.sum(residual**2)
        abs_losses[i] = np.sum(np.abs(residual))
    return sq_losses, abs_losses


def second_differences(values) -> np.ndarray:
    """The discrete second derivative of a sequence: ``diff(diff(values))``.

    Constant second differences mean the curve is a parabola -- smooth,
    with a single well-defined slope of the slope. Jumping second
    differences mean the curve bends only at particular points -- a
    piecewise-linear, kinked shape.
    """
    return np.diff(np.asarray(values, dtype=float), n=2)


# --------------------------------------------------------------------------
# 3. The normal equations: squared error's closed form
# --------------------------------------------------------------------------


def normal_equations(x, y):
    """Solve for (intercept, slope) directly from the normal equations.

    Squared error is smooth everywhere, so setting its derivative to zero
    gives a linear system to solve: ``(X^T X) beta = X^T y``. Absolute
    error has no derivative at a residual of zero, so no equivalent closed
    form exists for it -- fitting it requires an iterative solver instead
    (which is what ``QuantileRegressor`` runs).
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    design = np.column_stack([np.ones_like(x), x])
    beta = np.linalg.solve(design.T @ design, design.T @ y)
    return float(beta[0]), float(beta[1])


def fit_ols(x, y):
    """Fit ordinary least squares with scikit-learn; return (intercept, slope)."""
    model = LinearRegression().fit(np.asarray(x, dtype=float).reshape(-1, 1), y)
    return float(model.intercept_), float(model.coef_[0])


# --------------------------------------------------------------------------
# 4. Outlier sensitivity: squared error, Huber, and absolute error compared
# --------------------------------------------------------------------------


def fit_huber(x, y, epsilon: float = 1.35, max_iter: int = 500):
    """Fit scikit-learn's HuberRegressor; return (intercept, slope)."""
    model = HuberRegressor(epsilon=epsilon, max_iter=max_iter).fit(
        np.asarray(x, dtype=float).reshape(-1, 1), y
    )
    return float(model.intercept_), float(model.coef_[0])


def fit_quantile(x, y, quantile: float = 0.5):
    """Fit scikit-learn's QuantileRegressor at the median (absolute error's
    minimiser); return (intercept, slope). ``alpha=0`` disables the
    regulariser this estimator applies by default, so it measures plain
    absolute error.
    """
    model = QuantileRegressor(quantile=quantile, alpha=0.0, solver="highs").fit(
        np.asarray(x, dtype=float).reshape(-1, 1), y
    )
    return float(model.intercept_), float(model.coef_[0])


def outlier_shift(x, y, outlier_index: int | None = None, outlier_offset: float = 80.0):
    """Fit OLS, Huber and median regression before and after moving one
    point far off the line. Returns a dict with each estimator's slope
    before, slope after, and how far the slope moved.
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    if outlier_index is None:
        outlier_index = int(np.argmax(x))
    y_outlier = y.copy()
    y_outlier[outlier_index] = y_outlier[outlier_index] + outlier_offset

    result = {}
    for name, fitter in (("ols", fit_ols), ("huber", fit_huber), ("quantile", fit_quantile)):
        _b0, before = fitter(x, y)
        _a0, after = fitter(x, y_outlier)
        result[name] = {
            "before": round(before, 4),
            "after": round(after, 4),
            "movement": round(after - before, 4),
        }
    return result


def huber_epsilon_sweep(x, y, epsilons):
    """The Huber slope at each epsilon, holding the (outlier-contaminated)
    data fixed. Small epsilon leans on the absolute-error half of the
    loss and large epsilon leans on the squared-error half, converging to
    plain OLS as epsilon grows without bound.
    """
    rows = []
    for eps in epsilons:
        _intercept, slope = fit_huber(x, y, epsilon=eps)
        rows.append((float(eps), round(slope, 4)))
    return rows


# --------------------------------------------------------------------------
# 5. What each loss assumes: Gaussian errors, or something heavier-tailed
# --------------------------------------------------------------------------


def efficiency_under_noise(heavy_tailed: bool, replications: int = 500,
                            n: int = 150, true_slope: float = 3.0):
    """Fit OLS and Huber on many independent datasets with the same true
    line, and report the mean and standard deviation of each estimator's
    slope. Returns ``(ols_mean, ols_sd, huber_mean, huber_sd)``.

    Under Gaussian errors, Gauss-Markov says OLS is the best LINEAR
    UNBIASED estimator: among estimators that are linear in y and unbiased,
    OLS has the smallest variance. Under heavy-tailed errors that
    guarantee no longer implies OLS is the lowest-variance choice, and
    this function measures whether it still is.
    """
    ols_slopes = np.empty(replications)
    huber_slopes = np.empty(replications)
    for seed in range(replications):
        x, y = make_line_data(n=n, seed=seed, true_slope=true_slope, heavy_tailed=heavy_tailed)
        _i0, ols_slopes[seed] = fit_ols(x, y)
        _i1, huber_slopes[seed] = fit_huber(x, y)
    return (
        round(float(ols_slopes.mean()), 4),
        round(float(ols_slopes.std()), 4),
        round(float(huber_slopes.mean()), 4),
        round(float(huber_slopes.std()), 4),
    )

#!/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 loss_lib as L  # noqa: E402

VALUES = [2.0, 3.0, 5.0, 7.0, 100.0]


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


def main() -> None:
    print("Day 149 -- loss functions and least squares, measured")
    print("=" * 55)

    rule("1. What each loss minimises")
    best_sq = L.grid_minimize(VALUES, L.sse, 0.0, 110.0)
    best_abs = L.grid_minimize(VALUES, L.sae, 0.0, 110.0)
    print(f"  values: {VALUES}")
    print(f"  mean   = {np.mean(VALUES):.4f}   grid argmin of squared error = {best_sq:.5f}")
    print(f"  median = {np.median(VALUES):.4f}    grid argmin of absolute error = {best_abs:.5f}")

    rule("2. The shape of the loss landscape")
    x, y = L.make_line_data(n=40, seed=3)
    slopes = np.round(np.arange(2.0, 4.01, 0.1), 4)
    sq_losses, abs_losses = L.loss_landscape(x, y, intercept=5.0, slopes=slopes)
    print(f"  argmin, squared-error slope  : {slopes[int(np.argmin(sq_losses))]}")
    print(f"  argmin, absolute-error slope : {slopes[int(np.argmin(abs_losses))]}")
    print(f"  sd of squared-error's second differences  : {np.std(L.second_differences(sq_losses)):.6f} (constant -> parabola)")
    print(f"  sd of absolute-error's second differences : {np.std(L.second_differences(abs_losses)):.4f} (varies -> kinked)")

    rule("3. The normal equations, squared error's closed form")
    x3, y3 = L.make_line_data(n=300, seed=2)
    intercept_eq, slope_eq = L.normal_equations(x3, y3)
    intercept_sk, slope_sk = L.fit_ols(x3, y3)
    print(f"  normal equations : intercept={intercept_eq:.6f}  slope={slope_eq:.6f}")
    print(f"  LinearRegression : intercept={intercept_sk:.6f}  slope={slope_sk:.6f}")
    print(f"  difference       : intercept={abs(intercept_eq - intercept_sk):.2e}  slope={abs(slope_eq - slope_sk):.2e}")

    rule("4. Outlier sensitivity: one point, moved 80 units off the line")
    x4, y4 = L.make_line_data(n=60, seed=1, noise_sd=1.5)
    result = L.outlier_shift(x4, y4, outlier_offset=80.0)
    print("  estimator    before     after    movement")
    for name in ("ols", "huber", "quantile"):
        r = result[name]
        print(f"  {name:9s}  {r['before']:.4f}    {r['after']:.4f}    {r['movement']:+.4f}")
    ols_move = abs(result["ols"]["movement"])
    print(f"  OLS moved {ols_move / abs(result['huber']['movement']):.1f}x further than Huber, "
          f"{ols_move / abs(result['quantile']['movement']):.1f}x further than the median fit")

    rule("5. Huber's delta, sweeping from absolute-error-like to squared-error-like")
    y_outlier = np.asarray(y4, dtype=float).copy()
    y_outlier[int(np.argmax(x4))] += 80.0
    sweep = L.huber_epsilon_sweep(x4, y_outlier, [1.0, 1.35, 1.5, 2.0, 5.0, 20.0, 100.0])
    print("  epsilon      slope")
    for eps, slope in sweep:
        print(f"  {eps:7.2f}   {slope:.4f}")
    ols_intercept, ols_slope = L.fit_ols(x4, y_outlier)
    print(f"  OLS on the same outlier-contaminated data: slope={ols_slope:.4f} -- matches the large-epsilon end of the sweep")

    rule("6. What squared error assumes about the errors")
    gauss = L.efficiency_under_noise(heavy_tailed=False, replications=500)
    heavy = L.efficiency_under_noise(heavy_tailed=True, replications=500)
    print("  500 replications, n=150 rows each, true slope 3.0")
    print(f"  Gaussian errors     : OLS mean={gauss[0]:.4f} sd={gauss[1]:.4f}   Huber mean={gauss[2]:.4f} sd={gauss[3]:.4f}")
    print(f"    ratio sd(OLS)/sd(Huber) = {gauss[1] / gauss[3]:.4f}  (below 1: OLS is the more efficient choice)")
    print(f"  heavy-tailed errors : OLS mean={heavy[0]:.4f} sd={heavy[1]:.4f}   Huber mean={heavy[2]:.4f} sd={heavy[3]:.4f}")
    print(f"    ratio sd(OLS)/sd(Huber) = {heavy[1] / heavy[3]:.4f}  (above 1: Huber is now the more efficient choice)")
    print("  both estimators stay close to unbiased in both settings; only the spread changes")


if __name__ == "__main__":
    main()

"""Ten exercises in what choosing a loss function actually decides.

Read `../starter/00_brief.md` first. `loss_lib.py` is complete -- it is the
machinery, not the exercise. This file holds the reference solutions.

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 loss_lib as L

VALUES = [2.0, 3.0, 5.0, 7.0, 100.0]


def test_01_mean_minimizes_squared_error():
    best = L.grid_minimize(VALUES, L.sse, 0.0, 110.0)
    assert round(np.mean(VALUES), 4) == 23.4
    assert abs(best - 23.4) < 0.001


def test_01b_median_minimizes_absolute_error():
    best = L.grid_minimize(VALUES, L.sae, 0.0, 110.0)
    assert np.median(VALUES) == 5.0
    assert abs(best - 5.0) < 0.001
    # The mean is dragged toward the 100.0 outlier; the median ignores its
    # size entirely and only counts how many values sit on each side of it.
    assert np.mean(VALUES) > 4 * np.median(VALUES)


def test_02_squared_error_landscape_is_smooth_with_one_minimum():
    x, y = L.make_line_data(n=40, seed=3)
    slopes = np.round(np.arange(2.0, 4.01, 0.1), 4)
    sq_losses, _abs_losses = L.loss_landscape(x, y, intercept=5.0, slopes=slopes)
    assert slopes[int(np.argmin(sq_losses))] == 3.0
    # A parabola has a CONSTANT second difference; this is what "smooth,
    # one minimum" means numerically rather than just visually.
    second = L.second_differences(sq_losses)
    assert round(float(np.std(second)), 6) == 0.0


def test_02b_absolute_error_landscape_is_piecewise_linear_and_kinked():
    x, y = L.make_line_data(n=40, seed=3)
    slopes = np.round(np.arange(2.0, 4.01, 0.1), 4)
    _sq_losses, abs_losses = L.loss_landscape(x, y, intercept=5.0, slopes=slopes)
    assert slopes[int(np.argmin(abs_losses))] == 3.0
    # Absolute error's second difference is NOT constant: the slope of the
    # loss changes only at the slopes where a residual crosses zero, so the
    # curve is a sequence of straight segments meeting at kinks.
    second = L.second_differences(abs_losses)
    assert round(float(np.std(second)), 4) == 1.9366
    assert float(np.std(second)) > 0.0


def test_03_the_normal_equations_solve_squared_error_in_closed_form():
    x, y = L.make_line_data(n=300, seed=2)
    intercept_eq, slope_eq = L.normal_equations(x, y)
    intercept_sk, slope_sk = L.fit_ols(x, y)
    assert abs(intercept_eq - intercept_sk) < 1e-9
    assert abs(slope_eq - slope_sk) < 1e-9
    assert round(slope_eq, 4) == 2.9779
    assert round(intercept_eq, 4) == 4.9663


def test_04_ols_moves_far_when_a_single_point_becomes_an_outlier():
    x, y = L.make_line_data(n=60, seed=1, noise_sd=1.5)
    result = L.outlier_shift(x, y, outlier_offset=80.0)
    assert result["ols"] == {"before": 3.0465, "after": 3.801, "movement": 0.7545}


def test_04b_huber_and_median_regression_barely_move():
    x, y = L.make_line_data(n=60, seed=1, noise_sd=1.5)
    result = L.outlier_shift(x, y, outlier_offset=80.0)
    assert result["huber"] == {"before": 2.987, "after": 3.0308, "movement": 0.0437}
    assert result["quantile"] == {"before": 2.9961, "after": 3.0064, "movement": 0.0104}
    # OLS moved roughly 17x further than Huber and 72x further than the
    # median fit, from the identical single outlier.
    ols_move = abs(result["ols"]["movement"])
    assert round(ols_move / abs(result["huber"]["movement"]), 1) == 17.3
    assert round(ols_move / abs(result["quantile"]["movement"]), 1) == 72.5


def test_05_hubers_delta_interpolates_between_absolute_and_squared_error():
    x, y = L.make_line_data(n=60, seed=1, noise_sd=1.5)
    y_outlier = np.asarray(y, dtype=float).copy()
    y_outlier[int(np.argmax(x))] += 80.0
    sweep = L.huber_epsilon_sweep(x, y_outlier, [1.0, 1.35, 1.5, 2.0, 5.0, 20.0, 100.0])
    assert sweep == [
        (1.0, 3.0064),
        (1.35, 3.0308),
        (1.5, 3.0505),
        (2.0, 3.0906),
        (5.0, 3.1511),
        (20.0, 3.801),
        (100.0, 3.801),
    ]
    slopes = [s for _e, s in sweep]
    assert all(a <= b for a, b in zip(slopes, slopes[1:]))
    ols_intercept, ols_slope = L.fit_ols(x, y_outlier)
    assert sweep[-1][1] == round(ols_slope, 4)


def test_06_squared_error_is_the_most_efficient_choice_under_gaussian_errors():
    ols_mean, ols_sd, huber_mean, huber_sd = L.efficiency_under_noise(
        heavy_tailed=False, replications=500
    )
    assert (ols_mean, ols_sd, huber_mean, huber_sd) == (2.998, 0.056, 2.9977, 0.0588)
    assert abs(ols_mean - 3.0) < 0.01
    assert abs(huber_mean - 3.0) < 0.01
    # Under Gaussian errors OLS has the smaller spread -- Gauss-Markov's
    # promise that it is the BEST (lowest-variance) LINEAR UNBIASED
    # ESTIMATOR is not an abstraction here, it is this ratio being below 1.
    assert round(ols_sd / huber_sd, 4) == 0.9524
    assert ols_sd < huber_sd


def test_06b_but_not_under_heavy_tailed_errors():
    ols_mean, ols_sd, huber_mean, huber_sd = L.efficiency_under_noise(
        heavy_tailed=True, replications=500
    )
    assert (ols_mean, ols_sd, huber_mean, huber_sd) == (2.9967, 0.0589, 2.9984, 0.0422)
    assert abs(ols_mean - 3.0) < 0.01
    assert abs(huber_mean - 3.0) < 0.01
    # Both estimators are still roughly unbiased. But with fat-tailed
    # errors the ranking flips: Huber's spread is now the smaller one.
    assert round(ols_sd / huber_sd, 4) == 1.3957
    assert huber_sd < ols_sd

"""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 loss_lib as L


def test_sse_and_sae_are_zero_at_a_perfect_fit_and_positive_elsewhere():
    values = [1.0, 2.0, 3.0]
    assert L.sse(values, 2.0) == 2.0  # (1-2)^2 + (2-2)^2 + (3-2)^2
    assert L.sae(values, 2.0) == 2.0  # |1-2| + |2-2| + |3-2|
    assert L.sse(values, 100.0) > 0
    assert L.sae(values, 100.0) > 0


def test_make_line_data_is_deterministic_given_a_seed():
    x1, y1 = L.make_line_data(n=20, seed=42)
    x2, y2 = L.make_line_data(n=20, seed=42)
    assert np.array_equal(x1, x2)
    assert np.array_equal(y1, y2)
    x3, _y3 = L.make_line_data(n=20, seed=43)
    assert not np.array_equal(x1, x3)


def test_grid_minimize_finds_the_minimum_of_a_simple_quadratic():
    # sse(values, c) for a single value v is (v - c)^2, minimised at c = v.
    best = L.grid_minimize([7.0], L.sse, 0.0, 14.0, steps=14001)
    assert abs(best - 7.0) < 0.01


def test_normal_equations_recovers_a_known_line_exactly_when_there_is_no_noise():
    x = np.linspace(0.0, 10.0, 30)
    y = 5.0 + 3.0 * x  # no noise at all
    intercept, slope = L.normal_equations(x, y)
    assert abs(intercept - 5.0) < 1e-9
    assert abs(slope - 3.0) < 1e-9

lesson_id: D149
day: 149
kind: guided-build
languages:
  - python
  - bash
setup_commands:
  - cd labs/sections/machine-learning/day-149-loss-functions-and-least-squares
  - 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: 55
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 -> 14 passed.
  pytest starter -q -> 4 passed, 10 skipped (the four machinery checks in
  test_loss_lib.py are solved in both directories; the ten exercise stubs in
  starter/test_loss_claims.py are untouched). Everything ran through a real lab-local
  .venv created by the documented setup commands; scikit-learn pulled in scipy 1.18.1,
  joblib 1.5.3 and threadpoolctl 3.6.0 as its own dependencies, none of which this lab
  imports directly. The lab is fully offline after the pip install -- every dataset is
  a synthetic straight line generated on the spot from a seeded numpy.random.default_rng,
  nothing is downloaded, no dataset is bundled, and harness check 9 confirms no URL
  appears anywhere in starter/ or examples/ source. Section 7 of the harness copies
  examples/ into a mktemp -d scratch directory, confirms 14 passed, rewrites the
  exercise-4 OLS-movement assertion so its expected movement is 0.0 instead of 0.7545,
  confirms a non-zero exit naming the failing test, and removes the scratch directory.
  MEASURED PAIRS, all captured verbatim in expected-output/measured-values.txt. (1)
  WHAT EACH LOSS MINIMISES: on values [2.0, 3.0, 5.0, 7.0, 100.0], a fine grid search
  over sum((v-c)^2) lands at 23.39975 against the true mean of 23.4000, and the same
  search over sum(|v-c|) lands at 5.00005 against the true median of 5.0000 -- the
  single value of 100.0 drags the mean far from the other four points while the median
  ignores its size entirely. (2) LOSS-LANDSCAPE SHAPE: sweeping a candidate slope from
  2.0 to 4.0 on a 40-row dataset, squared error's second differences have standard
  deviation 0.000000 (a parabola: smooth, one minimum) while absolute error's have
  standard deviation 1.9366 (piecewise-linear, kinked at every point where a residual
  crosses zero); both landscapes are minimised at the identical slope of 3.0. (3) THE
  NORMAL EQUATIONS: solved directly with numpy.linalg.solve on a 300-row dataset, they
  give intercept=4.966335126423676 and slope=2.977853013652841, matching scikit-learn's
  LinearRegression (4.966335126423733, 2.977853013652829) to within 6e-14 on the
  intercept and 1e-14 on the slope -- squared error's smoothness is what makes this
  closed form possible; absolute error has no equivalent. (4) OUTLIER SENSITIVITY, the
  centrepiece: on a 60-row dataset, moving ONE point 80 units off the line moves
  LinearRegression's slope from 3.0465 to 3.8010 (+0.7545), HuberRegressor's from 2.9870
  to 3.0308 (+0.0437), and QuantileRegressor(quantile=0.5, alpha=0, solver='highs')'s
  from 2.9961 to 3.0064 (+0.0104) -- OLS moved 17.3x further than Huber and 72.5x
  further than the median fit. (5) HUBER'S DELTA SWEEP on the same outlier-contaminated
  data: epsilon 1.0 gives slope 3.0064, 1.35 (the scikit-learn default) gives 3.0308,
  1.5 gives 3.0505, 2.0 gives 3.0906, 5.0 gives 3.1511, and both 20.0 and 100.0 give
  3.8010, exactly matching plain OLS on the same data -- the sweep is monotonic and
  converges exactly to least squares as epsilon grows. (6) WHAT SQUARED ERROR ASSUMES:
  over 500 replications of 150-row datasets with a true slope of 3.0, Gaussian errors
  give OLS mean=2.998 sd=0.0560 against Huber mean=2.9977 sd=0.0588 (ratio 0.9524,
  OLS the tighter estimator), while Student's-t (df=3) heavy-tailed errors of similar
  central spread give OLS mean=2.9967 sd=0.0589 against Huber mean=2.9984 sd=0.0422
  (ratio 1.3957, Huber now the tighter estimator) -- both estimators stay within 0.01
  of the true slope in both settings, so only the precision ranking flips, not the bias.
  Harness check 8 re-confirms the outlier-sensitivity ranking at five dataset seeds, the
  normal equations at three further sample sizes, and the Gauss-Markov ranking at a
  shorter replication count (150), so no directional claim in this lab rests on the
  single seed or replication count quoted in the lesson. TWO HONESTY CALLS. FIRST:
  HuberRegressor and QuantileRegressor were both verified to exist and converge in
  scikit-learn 1.9.0, without a warning, on every dataset used in this lab, before
  anything was built on them -- neither was assumed from memory or from an older
  version's documentation. QuantileRegressor is fit with alpha=0.0 and solver='highs'
  explicitly, because its default alpha applies an L1 penalty that this lab's plain
  absolute-error claim does not want mixed in, and naming the solver avoids depending
  on whichever default happens to be current. SECOND: Huber's mean slope under both
  error settings (2.9977 and 2.9984) sits slightly further from the true value of 3.0
  than OLS's Gaussian-error mean (2.998) does, which is sampling noise in an average of
  500 draws rather than evidence that Huber is biased -- harness check 8 confirms both
  estimators' means stay within 0.01 of the truth at a different replication count too,
  and no claim in this lab treats one estimator's mean as more "correct" than the
  other's; the whole point of exercises 6 and 6b is that both are close to unbiased and
  only their spread differs, and which spread is smaller depends on the errors.

# 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 sampled figure in this lab — the outlier-shift movements, the Huber
epsilon sweep on real data, and both efficiency-under-noise comparisons —
comes from `numpy.random.default_rng`, and NumPy's own documentation is
explicit that `Generator` makes no promise of stream compatibility between
versions. A different NumPy can legitimately produce a different stream
from the same seed, and every sampled figure would move.

`HuberRegressor` and `QuantileRegressor` are also estimator internals: the
exact coefficients from an iterative solver can shift in the last few
decimal places across scikit-learn releases even with the seed held fixed,
because the solver's stopping tolerance and default number of iterations
are part of the library, not the seed.

What does not depend on the pins: the mean minimises squared error and the
median minimises absolute error, which is arithmetic; the normal equations
matching `LinearRegression` to many decimal places, which is closed-form
algebra; the constant-versus-varying second differences that distinguish a
smooth loss from a kinked one; and the direction of every comparison —
least squares moves further than Huber and further still than the median
fit when a single point becomes an outlier, and Gaussian errors favour
least squares while heavy-tailed errors favour Huber.
`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: every
dataset in this lab is generated on the spot from a seeded generator, and
nothing is downloaded.

## Free and open-source status

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

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

# Day 149 lab brief — Loss Functions and Least Squares

A loss function is not part of the mathematics of a straight line. It is a
choice you make, and this lab measures what that choice actually decides.

## The claim you are here to measure

> Squared error and absolute error agree on the easy cases and disagree
> sharply on the hard ones, and the disagreement is exactly outlier
> sensitivity.

Exercise 1 starts with five numbers: `2.0, 3.0, 5.0, 7.0, 100.0`. A grid
search over every candidate value confirms two textbook facts as
measurements rather than assertions:

| loss | minimiser | value here |
| --- | --- | --- |
| squared error, `sum((v - c)^2)` | the **mean** | 23.4000 |
| absolute error, `sum(|v - c|)` | the **median** | 5.0000 |

The mean is dragged nearly to 23.4 by the single value of 100.0. The
median does not move — it only counts how many points sit on each side of
it, never how far away they are. That single fact is the whole reason a
squared-error fit and an absolute-error fit disagree about outliers.

## The shape of the two landscapes

Exercise 2 sweeps a candidate slope from 2.0 to 4.0 against the same
sixty-row dataset and totals both losses at each value:

| loss | second differences | what that means |
| --- | --- | --- |
| squared error | constant (sd = 0.000000) | a parabola: smooth, one minimum |
| absolute error | varies (sd = 1.9366) | piecewise-linear: kinked at every point where a residual crosses zero |

Constant curvature is why squared error has a closed-form solution.
Exercise 3 solves that closed form directly — the **normal equations** —
and confirms it matches `LinearRegression` to thirteen decimal places.
Absolute error has no equivalent closed form, because its derivative does
not exist at a residual of exactly zero; fitting it needs an iterative
solver instead.

## The centrepiece: one outlier, three losses

Exercises 4 and 4b move a single point 80 units off the line and re-fit
three estimators on identical data:

| estimator | loss it minimises | slope before | slope after | movement |
| --- | --- | --- | --- | --- |
| `LinearRegression` | squared error | 3.0465 | 3.8010 | **+0.7545** |
| `HuberRegressor` | Huber (blended) | 2.9870 | 3.0308 | +0.0437 |
| `QuantileRegressor(0.5)` | absolute error | 2.9961 | 3.0064 | +0.0104 |

One point, and least squares moves seventeen times further than Huber and
seventy-two times further than the median fit. Squaring the residual
means a residual of 80 contributes 6,400 to the total loss instead of 80
— the optimizer will trade a great deal of fit elsewhere to shrink that
one huge term, and "elsewhere" is every other point on the line.

Exercise 5 sweeps Huber's `epsilon` parameter on that same contaminated
data, from 1.0 up to 100.0. The slope climbs smoothly from
absolute-error-like behaviour toward the least-squares answer, and at
large epsilon it lands exactly on the OLS slope — Huber's blend of the two
losses is genuinely continuous, not a discrete switch.

## What squared error is betting on

Exercises 6 and 6b are the payoff. Fit OLS and Huber on 500 independently
generated datasets, twice — once with Gaussian errors, once with
heavy-tailed (Student's t, 3 degrees of freedom) errors of similar central
spread:

| errors | OLS spread (sd) | Huber spread (sd) | ratio OLS/Huber |
| --- | --- | --- | --- |
| Gaussian | 0.0560 | 0.0588 | 0.9524 — OLS tighter |
| heavy-tailed | 0.0589 | 0.0422 | 1.3957 — Huber tighter |

Both estimators stay close to unbiased in both settings. What changes is
which one is more *precise*. Gauss-Markov's theorem says ordinary least
squares is the **best linear unbiased estimator** under its assumptions —
every one of those four words is load-bearing, and "best" specifically
means lowest variance *among linear, unbiased estimators*, not lowest
variance full stop, and not under any distribution of errors whatsoever.
Change the error distribution and the ranking measurably flips.

## 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_loss_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 `loss_lib.py`, `test_loss_lib.py` and `test_loss_claims.py`; pytest
aborts on the module-name collision. Run them separately, always.

"""Loss functions, measured: what choosing squared error over absolute error
actually decides, and what each is implicitly betting on.

A loss function is a choice, not a law of nature. This module measures four
consequences of that choice: the minimiser of squared error is the mean and
the minimiser of absolute error is the median; the squared-error landscape is
smooth with one minimum while the absolute-error landscape is piecewise
linear and kinked; the normal equations solve squared error in closed form
because that smoothness gives a zero derivative one can solve for directly,
while absolute error has no such closed form; and swapping which loss you
minimise changes how far a single outlier can move your line, and which loss
wins depends on what the errors actually look like -- Gaussian, or
heavy-tailed.

Everything here is deterministic given a seed.
"""

from __future__ import annotations

import numpy as np

from sklearn.linear_model import HuberRegressor, LinearRegression, QuantileRegressor


# --------------------------------------------------------------------------
# 1. What each loss minimises
# --------------------------------------------------------------------------


def sse(values, candidate: float) -> float:
    """Sum of squared error between a scalar candidate and every value."""
    values = np.asarray(values, dtype=float)
    return float(np.sum((values - candidate) ** 2))


def sae(values, candidate: float) -> float:
    """Sum of absolute error between a scalar candidate and every value."""
    values = np.asarray(values, dtype=float)
    return float(np.sum(np.abs(values - candidate)))


def grid_minimize(values, loss_fn, lo: float, hi: float, steps: int = 200_001) -> float:
    """Find the candidate in a fine grid that minimises the given loss.

    A brute-force numerical stand-in for calculus: since it searches a
    finite grid it lands within ``(hi - lo) / (steps - 1)`` of the true
    minimiser, not exactly on it.
    """
    grid = np.linspace(lo, hi, steps)
    losses = np.array([loss_fn(values, c) for c in grid])
    return float(grid[int(np.argmin(losses))])


# --------------------------------------------------------------------------
# 2. The shape of the loss landscape: smooth, or kinked
# --------------------------------------------------------------------------


def make_line_data(n: int = 150, seed: int = 0, true_intercept: float = 5.0,
                    true_slope: float = 3.0, noise_sd: float = 2.0,
                    heavy_tailed: bool = False, heavy_df: int = 3,
                    heavy_scale: float = 1.2):
    """A simple one-predictor dataset with a known true line.

    ``heavy_tailed=True`` swaps the Gaussian error for a scaled Student's t
    with ``heavy_df`` degrees of freedom, which has the same rough central
    spread but far fatter tails -- the construction used to measure what
    squared error implicitly assumes about the errors.
    """
    rng = np.random.default_rng(seed)
    x = rng.uniform(0.0, 10.0, n)
    if heavy_tailed:
        errors = rng.standard_t(df=heavy_df, size=n) * heavy_scale
    else:
        errors = rng.normal(0.0, noise_sd, n)
    y = true_intercept + true_slope * x + errors
    return x, y


def loss_landscape(x, y, intercept: float, slopes):
    """Total squared error and total absolute error at each candidate slope,
    with the intercept held fixed. Returns ``(sq_losses, abs_losses)``.
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    slopes = np.asarray(slopes, dtype=float)
    sq_losses = np.empty_like(slopes)
    abs_losses = np.empty_like(slopes)
    for i, m in enumerate(slopes):
        residual = y - (intercept + m * x)
        sq_losses[i] = np.sum(residual**2)
        abs_losses[i] = np.sum(np.abs(residual))
    return sq_losses, abs_losses


def second_differences(values) -> np.ndarray:
    """The discrete second derivative of a sequence: ``diff(diff(values))``.

    Constant second differences mean the curve is a parabola -- smooth,
    with a single well-defined slope of the slope. Jumping second
    differences mean the curve bends only at particular points -- a
    piecewise-linear, kinked shape.
    """
    return np.diff(np.asarray(values, dtype=float), n=2)


# --------------------------------------------------------------------------
# 3. The normal equations: squared error's closed form
# --------------------------------------------------------------------------


def normal_equations(x, y):
    """Solve for (intercept, slope) directly from the normal equations.

    Squared error is smooth everywhere, so setting its derivative to zero
    gives a linear system to solve: ``(X^T X) beta = X^T y``. Absolute
    error has no derivative at a residual of zero, so no equivalent closed
    form exists for it -- fitting it requires an iterative solver instead
    (which is what ``QuantileRegressor`` runs).
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    design = np.column_stack([np.ones_like(x), x])
    beta = np.linalg.solve(design.T @ design, design.T @ y)
    return float(beta[0]), float(beta[1])


def fit_ols(x, y):
    """Fit ordinary least squares with scikit-learn; return (intercept, slope)."""
    model = LinearRegression().fit(np.asarray(x, dtype=float).reshape(-1, 1), y)
    return float(model.intercept_), float(model.coef_[0])


# --------------------------------------------------------------------------
# 4. Outlier sensitivity: squared error, Huber, and absolute error compared
# --------------------------------------------------------------------------


def fit_huber(x, y, epsilon: float = 1.35, max_iter: int = 500):
    """Fit scikit-learn's HuberRegressor; return (intercept, slope)."""
    model = HuberRegressor(epsilon=epsilon, max_iter=max_iter).fit(
        np.asarray(x, dtype=float).reshape(-1, 1), y
    )
    return float(model.intercept_), float(model.coef_[0])


def fit_quantile(x, y, quantile: float = 0.5):
    """Fit scikit-learn's QuantileRegressor at the median (absolute error's
    minimiser); return (intercept, slope). ``alpha=0`` disables the
    regulariser this estimator applies by default, so it measures plain
    absolute error.
    """
    model = QuantileRegressor(quantile=quantile, alpha=0.0, solver="highs").fit(
        np.asarray(x, dtype=float).reshape(-1, 1), y
    )
    return float(model.intercept_), float(model.coef_[0])


def outlier_shift(x, y, outlier_index: int | None = None, outlier_offset: float = 80.0):
    """Fit OLS, Huber and median regression before and after moving one
    point far off the line. Returns a dict with each estimator's slope
    before, slope after, and how far the slope moved.
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    if outlier_index is None:
        outlier_index = int(np.argmax(x))
    y_outlier = y.copy()
    y_outlier[outlier_index] = y_outlier[outlier_index] + outlier_offset

    result = {}
    for name, fitter in (("ols", fit_ols), ("huber", fit_huber), ("quantile", fit_quantile)):
        _b0, before = fitter(x, y)
        _a0, after = fitter(x, y_outlier)
        result[name] = {
            "before": round(before, 4),
            "after": round(after, 4),
            "movement": round(after - before, 4),
        }
    return result


def huber_epsilon_sweep(x, y, epsilons):
    """The Huber slope at each epsilon, holding the (outlier-contaminated)
    data fixed. Small epsilon leans on the absolute-error half of the
    loss and large epsilon leans on the squared-error half, converging to
    plain OLS as epsilon grows without bound.
    """
    rows = []
    for eps in epsilons:
        _intercept, slope = fit_huber(x, y, epsilon=eps)
        rows.append((float(eps), round(slope, 4)))
    return rows


# --------------------------------------------------------------------------
# 5. What each loss assumes: Gaussian errors, or something heavier-tailed
# --------------------------------------------------------------------------


def efficiency_under_noise(heavy_tailed: bool, replications: int = 500,
                            n: int = 150, true_slope: float = 3.0):
    """Fit OLS and Huber on many independent datasets with the same true
    line, and report the mean and standard deviation of each estimator's
    slope. Returns ``(ols_mean, ols_sd, huber_mean, huber_sd)``.

    Under Gaussian errors, Gauss-Markov says OLS is the best LINEAR
    UNBIASED estimator: among estimators that are linear in y and unbiased,
    OLS has the smallest variance. Under heavy-tailed errors that
    guarantee no longer implies OLS is the lowest-variance choice, and
    this function measures whether it still is.
    """
    ols_slopes = np.empty(replications)
    huber_slopes = np.empty(replications)
    for seed in range(replications):
        x, y = make_line_data(n=n, seed=seed, true_slope=true_slope, heavy_tailed=heavy_tailed)
        _i0, ols_slopes[seed] = fit_ols(x, y)
        _i1, huber_slopes[seed] = fit_huber(x, y)
    return (
        round(float(ols_slopes.mean()), 4),
        round(float(ols_slopes.std()), 4),
        round(float(huber_slopes.mean()), 4),
        round(float(huber_slopes.std()), 4),
    )

"""Ten exercises in what choosing a loss function actually decides.

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.
`loss_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 loss_lib as L  # noqa: F401  (you will need it)

VALUES = [2.0, 3.0, 5.0, 7.0, 100.0]


def test_01_mean_minimizes_squared_error():
    pytest.skip(
        "Call L.grid_minimize(VALUES, L.sse, 0.0, 110.0) and assert it is "
        "within 0.001 of np.mean(VALUES), which rounds to 23.4. The mean is "
        "not chosen for squared error by convention -- it is what a grid "
        "search over the loss actually finds."
    )


def test_01b_median_minimizes_absolute_error():
    pytest.skip(
        "Call L.grid_minimize(VALUES, L.sae, 0.0, 110.0) and assert it is "
        "within 0.001 of np.median(VALUES), which is 5.0. Then assert "
        "np.mean(VALUES) is more than four times np.median(VALUES): the "
        "single value of 100.0 drags the mean far from the other four "
        "points, while the median only counts how many values sit on each "
        "side of it and ignores their size entirely."
    )


def test_02_squared_error_landscape_is_smooth_with_one_minimum():
    pytest.skip(
        "Build x, y = L.make_line_data(n=40, seed=3). Sweep slopes = "
        "np.round(np.arange(2.0, 4.01, 0.1), 4) with the intercept fixed "
        "at 5.0 using L.loss_landscape, and assert the squared-error "
        "minimiser is exactly slope 3.0. Then assert "
        "round(np.std(L.second_differences(sq_losses)), 6) == 0.0 -- a "
        "CONSTANT second difference is the numerical signature of a "
        "parabola: smooth, with one minimum."
    )


def test_02b_absolute_error_landscape_is_piecewise_linear_and_kinked():
    pytest.skip(
        "Using the same x, y and slopes as test_02, assert the "
        "absolute-error minimiser is also exactly slope 3.0, but that "
        "round(np.std(L.second_differences(abs_losses)), 4) == 1.9366 -- "
        "NOT constant. Absolute error only bends where a residual crosses "
        "zero, so its landscape is a sequence of straight segments meeting "
        "at kinks rather than one smooth curve."
    )


def test_03_the_normal_equations_solve_squared_error_in_closed_form():
    pytest.skip(
        "Build x, y = L.make_line_data(n=300, seed=2). Solve with "
        "L.normal_equations(x, y) and separately with L.fit_ols(x, y), and "
        "assert the two intercepts agree to within 1e-9 and the two slopes "
        "agree to within 1e-9. Then assert the normal-equations slope "
        "rounds to 2.9779 and the intercept to 4.9663. Squared error is "
        "smooth everywhere, so setting its derivative to zero gives a "
        "linear system with an exact solution -- absolute error has no "
        "equivalent closed form, because it has no derivative at a "
        "residual of exactly zero."
    )


def test_04_ols_moves_far_when_a_single_point_becomes_an_outlier():
    pytest.skip(
        "Build x, y = L.make_line_data(n=60, seed=1, noise_sd=1.5) and call "
        "L.outlier_shift(x, y, outlier_offset=80.0). Assert result['ols'] "
        "equals {'before': 3.0465, 'after': 3.801, 'movement': 0.7545}. "
        "Moving ONE point 80 units off the line moves the least-squares "
        "slope by three quarters of a unit."
    )


def test_04b_huber_and_median_regression_barely_move():
    pytest.skip(
        "Using the same result from test_04, assert result['huber'] equals "
        "{'before': 2.987, 'after': 3.0308, 'movement': 0.0437} and "
        "result['quantile'] equals {'before': 2.9961, 'after': 3.0064, "
        "'movement': 0.0104}. Then assert OLS's movement, divided by "
        "Huber's, rounds to 17.3, and divided by the median fit's, rounds "
        "to 72.5. Same data, same outlier, same single point moved -- only "
        "the loss changed."
    )


def test_05_hubers_delta_interpolates_between_absolute_and_squared_error():
    pytest.skip(
        "Reuse x, y from test_04, copy y, add 80.0 to the row at "
        "np.argmax(x), and call L.huber_epsilon_sweep on the contaminated "
        "data with epsilons [1.0, 1.35, 1.5, 2.0, 5.0, 20.0, 100.0]. Assert "
        "the result equals [(1.0, 3.0064), (1.35, 3.0308), (1.5, 3.0505), "
        "(2.0, 3.0906), (5.0, 3.1511), (20.0, 3.801), (100.0, 3.801)]. "
        "Assert the slopes are non-decreasing in epsilon, and that the "
        "final slope equals L.fit_ols on the same contaminated data, "
        "rounded to 4 places. Small epsilon leans on absolute error; large "
        "epsilon converges to plain least squares."
    )


def test_06_squared_error_is_the_most_efficient_choice_under_gaussian_errors():
    pytest.skip(
        "Call L.efficiency_under_noise(heavy_tailed=False, "
        "replications=500) and assert it equals (2.998, 0.056, 2.9977, "
        "0.0588) -- (ols_mean, ols_sd, huber_mean, huber_sd). Assert both "
        "means are within 0.01 of the true slope of 3.0 (both estimators "
        "are roughly unbiased), then assert round(ols_sd / huber_sd, 4) == "
        "0.9524 and that ols_sd < huber_sd. Under Gaussian errors, OLS has "
        "the smaller spread -- this ratio being below 1 is what Gauss-"
        "Markov's promise of the BEST (lowest-variance) LINEAR UNBIASED "
        "ESTIMATOR looks like when you measure it."
    )


def test_06b_but_not_under_heavy_tailed_errors():
    pytest.skip(
        "Call L.efficiency_under_noise(heavy_tailed=True, "
        "replications=500) and assert it equals (2.9967, 0.0589, 2.9984, "
        "0.0422). Assert both means are still within 0.01 of 3.0, then "
        "assert round(ols_sd / huber_sd, 4) == 1.3957 and that huber_sd < "
        "ols_sd -- the ranking has FLIPPED. Gauss-Markov's guarantee is "
        "conditional on the errors; change what the errors look like and "
        "the best linear unbiased estimator is no longer the most precise "
        "one available."
    )

"""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 loss_lib as L


def test_sse_and_sae_are_zero_at_a_perfect_fit_and_positive_elsewhere():
    values = [1.0, 2.0, 3.0]
    assert L.sse(values, 2.0) == 2.0  # (1-2)^2 + (2-2)^2 + (3-2)^2
    assert L.sae(values, 2.0) == 2.0  # |1-2| + |2-2| + |3-2|
    assert L.sse(values, 100.0) > 0
    assert L.sae(values, 100.0) > 0


def test_make_line_data_is_deterministic_given_a_seed():
    x1, y1 = L.make_line_data(n=20, seed=42)
    x2, y2 = L.make_line_data(n=20, seed=42)
    assert np.array_equal(x1, x2)
    assert np.array_equal(y1, y2)
    x3, _y3 = L.make_line_data(n=20, seed=43)
    assert not np.array_equal(x1, x3)


def test_grid_minimize_finds_the_minimum_of_a_simple_quadratic():
    # sse(values, c) for a single value v is (v - c)^2, minimised at c = v.
    best = L.grid_minimize([7.0], L.sse, 0.0, 14.0, steps=14001)
    assert abs(best - 7.0) < 0.01


def test_normal_equations_recovers_a_known_line_exactly_when_there_is_no_noise():
    x = np.linspace(0.0, 10.0, 30)
    y = 5.0 + 3.0 * x  # no noise at all
    intercept, slope = L.normal_equations(x, y)
    assert abs(intercept - 5.0) < 1e-9
    assert abs(slope - 3.0) < 1e-9

#!/usr/bin/env bash
# Day 149 lab harness: "Loss Functions and Least Squares"
#
# Prints "N checks, M failure(s)" and exits 0 only when M is zero.
set -u

LAB_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)"
cd "$LAB_DIR"

PYTHON="${PYTHON:-.venv/bin/python3}"
PYTEST="${PYTEST:-.venv/bin/pytest}"

# Clear caches at the START so the final cleanliness check measures what
# THIS run left behind, not what a previous `pytest starter -q` left.
find . -path ./.venv -prune -o -type d -name '__pycache__' -exec rm -rf -- {} + 2>/dev/null
rm -rf .pytest_cache

CHECKS=0
FAILURES=0

ok() {
  CHECKS=$((CHECKS + 1))
  echo "  ok: $1"
}

fail() {
  CHECKS=$((CHECKS + 1))
  FAILURES=$((FAILURES + 1))
  echo "  FAIL: $1"
}

if [ ! -x "$PYTHON" ]; then
  echo "No lab .venv found at $PYTHON."
  echo "Run: python3 -m venv .venv && .venv/bin/pip install -r requirements/requirements.txt"
  exit 2
fi

echo "1. Installed versions match requirements/requirements.txt"
VERSION_CHECK=$("$PYTHON" - <<'PYEOF'
import numpy, sklearn, pytest
print("numpy", numpy.__version__)
print("scikit-learn", sklearn.__version__)
print("pytest", pytest.__version__)
PYEOF
)
echo "$VERSION_CHECK" | sed 's/^/    /'
while read -r pkg pin; do
  pin_version="${pin#*==}"
  installed=$(echo "$VERSION_CHECK" | awk -v p="$pkg" '$1==p {print $2}')
  if [ "$installed" = "$pin_version" ]; then
    ok "$pkg $installed matches the pin"
  else
    fail "$pkg installed=$installed pinned=$pin_version"
  fi
done < <(sed 's/==/ ==/' requirements/requirements.txt)

echo ""
echo "2. Every published claim, reproduced directly (no pytest involved)"
DIRECT_CHECK=$("$PYTHON" - <<'PYEOF'
import sys
sys.path.insert(0, "examples")

import numpy as np

import loss_lib as L

errors = []
VALUES = [2.0, 3.0, 5.0, 7.0, 100.0]


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


# 1. Mean and median
best_sq = L.grid_minimize(VALUES, L.sse, 0.0, 110.0)
best_abs = L.grid_minimize(VALUES, L.sae, 0.0, 110.0)
expect("mean", round(float(np.mean(VALUES)), 4), 23.4)
if abs(best_sq - 23.4) >= 0.001:
    errors.append(f"grid argmin of squared error {best_sq} not close to the mean")
expect("median", float(np.median(VALUES)), 5.0)
if abs(best_abs - 5.0) >= 0.001:
    errors.append(f"grid argmin of absolute error {best_abs} not close to the median")

# 2. Landscape shape
x, y = L.make_line_data(n=40, seed=3)
slopes = np.round(np.arange(2.0, 4.01, 0.1), 4)
sq_losses, abs_losses = L.loss_landscape(x, y, intercept=5.0, slopes=slopes)
expect("squared-error argmin slope", float(slopes[int(np.argmin(sq_losses))]), 3.0)
expect("absolute-error argmin slope", float(slopes[int(np.argmin(abs_losses))]), 3.0)
expect("sd of squared-error second differences", round(float(np.std(L.second_differences(sq_losses))), 6), 0.0)
expect("sd of absolute-error second differences", round(float(np.std(L.second_differences(abs_losses))), 4), 1.9366)

# 3. Normal equations
x3, y3 = L.make_line_data(n=300, seed=2)
intercept_eq, slope_eq = L.normal_equations(x3, y3)
intercept_sk, slope_sk = L.fit_ols(x3, y3)
if abs(intercept_eq - intercept_sk) >= 1e-9:
    errors.append(f"normal-equations intercept {intercept_eq} vs sklearn {intercept_sk}")
if abs(slope_eq - slope_sk) >= 1e-9:
    errors.append(f"normal-equations slope {slope_eq} vs sklearn {slope_sk}")
expect("normal-equations slope", round(slope_eq, 4), 2.9779)
expect("normal-equations intercept", round(intercept_eq, 4), 4.9663)

# 4. Outlier shift
x4, y4 = L.make_line_data(n=60, seed=1, noise_sd=1.5)
result = L.outlier_shift(x4, y4, outlier_offset=80.0)
expect("ols outlier shift", result["ols"], {"before": 3.0465, "after": 3.801, "movement": 0.7545})
expect("huber outlier shift", result["huber"], {"before": 2.987, "after": 3.0308, "movement": 0.0437})
expect("quantile outlier shift", result["quantile"], {"before": 2.9961, "after": 3.0064, "movement": 0.0104})
ols_move = abs(result["ols"]["movement"])
expect("ols/huber movement ratio", round(ols_move / abs(result["huber"]["movement"]), 1), 17.3)
expect("ols/quantile movement ratio", round(ols_move / abs(result["quantile"]["movement"]), 1), 72.5)

# 5. Huber epsilon sweep
y_outlier = np.asarray(y4, dtype=float).copy()
y_outlier[int(np.argmax(x4))] += 80.0
sweep = L.huber_epsilon_sweep(x4, y_outlier, [1.0, 1.35, 1.5, 2.0, 5.0, 20.0, 100.0])
expect(
    "huber epsilon sweep",
    sweep,
    [(1.0, 3.0064), (1.35, 3.0308), (1.5, 3.0505), (2.0, 3.0906), (5.0, 3.1511), (20.0, 3.801), (100.0, 3.801)],
)
sweep_slopes = [s for _e, s in sweep]
if not all(a <= b for a, b in zip(sweep_slopes, sweep_slopes[1:])):
    errors.append("huber epsilon sweep was not non-decreasing")
_i, ols_slope_outlier = L.fit_ols(x4, y_outlier)
expect("sweep converges to OLS", sweep[-1][1], round(ols_slope_outlier, 4))

# 6. Efficiency under noise
gauss = L.efficiency_under_noise(heavy_tailed=False, replications=500)
heavy = L.efficiency_under_noise(heavy_tailed=True, replications=500)
expect("gaussian efficiency", gauss, (2.998, 0.056, 2.9977, 0.0588))
expect("heavy-tailed efficiency", heavy, (2.9967, 0.0589, 2.9984, 0.0422))
if not (gauss[1] < gauss[3]):
    errors.append("OLS was not the tighter estimator under Gaussian errors")
if not (heavy[3] < heavy[1]):
    errors.append("Huber was not the tighter estimator under heavy-tailed errors")
expect("ratio under gaussian errors", round(gauss[1] / gauss[3], 4), 0.9524)
expect("ratio under heavy-tailed errors", round(heavy[1] / heavy[3], 4), 1.3957)

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-6 reproduced directly against loss_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 "^14 passed"; then
  ok "pytest examples -q -> 14 passed"
else
  fail "pytest examples -q did not report 14 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, 10 skipped"; then
  ok "pytest starter -q -> 4 passed, 10 skipped (the machinery checks pass; the ten exercises are stubs)"
else
  fail "pytest starter -q did not report 4 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}/d149-scratch.XXXXXX")
cp examples/*.py "$SCRATCH"/
SCRATCH_OUT=$("$PYTEST" "$SCRATCH" -q 2>&1)
if echo "$SCRATCH_OUT" | tail -1 | grep -qE "^14 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_loss_claims.py" <<'PYEOF'
import sys
path = sys.argv[1]
text = open(path).read()
needle = "assert result[\"ols\"] == {\"before\": 3.0465, \"after\": 3.801, \"movement\": 0.7545}"
replacement = "assert result[\"ols\"] == {\"before\": 3.0465, \"after\": 3.801, \"movement\": 0.0}"
assert needle in text, "could not find the assertion to break"
open(path, "w").write(text.replace(needle, replacement, 1))
PYEOF
BROKEN_OUT=$("$PYTEST" "$SCRATCH" -q 2>&1)
BROKEN_STATUS=$?
if [ "$BROKEN_STATUS" -ne 0 ] && echo "$BROKEN_OUT" | grep -q "test_04_ols_moves_far_when_a_single_point_becomes_an_outlier"; then
  ok "breaking exercise 4'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 seeds"
DIRECTION=$("$PYTHON" - <<'PYEOF'
import sys
sys.path.insert(0, "examples")
import numpy as np
import loss_lib as L

problems = []

# Outlier sensitivity is not a property of one dataset seed.
for seed in (1, 2, 3, 4, 5):
    x, y = L.make_line_data(n=60, seed=seed, noise_sd=1.5)
    result = L.outlier_shift(x, y, outlier_offset=80.0)
    ols_move = abs(result["ols"]["movement"])
    huber_move = abs(result["huber"]["movement"])
    quantile_move = abs(result["quantile"]["movement"])
    if not (ols_move > huber_move and ols_move > quantile_move):
        problems.append(f"seed {seed}: OLS did not move furthest (ols={ols_move}, huber={huber_move}, quantile={quantile_move})")

# The normal equations match LinearRegression at other sample sizes too.
for n, seed in ((50, 10), (500, 11), (1000, 12)):
    x, y = L.make_line_data(n=n, seed=seed)
    i_eq, s_eq = L.normal_equations(x, y)
    i_sk, s_sk = L.fit_ols(x, y)
    if abs(i_eq - i_sk) >= 1e-8 or abs(s_eq - s_sk) >= 1e-8:
        problems.append(f"n={n} seed={seed}: normal equations diverged from sklearn")

# Gauss-Markov's ranking is not a property of one replication count.
gauss_short = L.efficiency_under_noise(heavy_tailed=False, replications=150)
heavy_short = L.efficiency_under_noise(heavy_tailed=True, replications=150)
if not (gauss_short[1] < gauss_short[3]):
    problems.append("OLS was not tighter than Huber under Gaussian errors at 150 replications")
if not (heavy_short[3] < heavy_short[1]):
    problems.append("Huber was not tighter than OLS under heavy-tailed errors at 150 replications")

if problems:
    for p in problems:
        print("ERROR:", p)
else:
    print("every direction held")
PYEOF
)
if [ "$DIRECTION" = "every direction held" ]; then
  ok "outlier sensitivity, the normal equations and the Gauss-Markov ranking hold at seeds the lesson does not quote"
else
  fail "a direction failed beyond the quoted seed"
  echo "$DIRECTION" | sed 's/^/    /'
fi

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

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